Hi all. Thanks a lot for developing and maintaining the pipeline! I would like to perform the alignment using the GATK.GRCh38 reference genome to maintain consistency with the subsequent variant calling step and improving alignment accuracy for variant calling (as it include extra decoy and ALT contigs). Does it make sense to use this specific genome for RNA-seq alignment? If so, where can I find a suitable GTF compatible with the GATK.GRCh38 reference, considering the added contigs? I've checked the GATK resource bundle, but it doesn't seem to include any annotation file. I'm concerned that a standard GTF from sources like GENCODE might not fully account for the extra contigs, potentially causing issues. Thanks in advance!

You certainly can use the GATK genome for RNAseq, but I'm not sure it is optimal, because extra haplotype sequences could cause unnecessary multimapping. One issue you are likely to run into is that the GATK genome uses chromosome names with the

chr

prefix, whereas most other versions (Ensembl, GENCODE, NCBI) use just the straight number/letter (

1, 2, 3, ..., 22, X, Y, MT

). But for the primary contigs there are "translation" tables, if you need, and the vast majority of genes are on the "canonical" chromosomes, so just removing/adding

chr

will do the trick. (If you have WES, all of your variant calls will be on canonical chromosomes, and if you use use the "standard" intervals file for WGS, the same is true there. The GATK genome literally has thousands of contigs, but a lot of those are decoys from other species, designed to attracted "garbage" reads and then be ignored. And even the "random" contigs from the human genome don't have a lot of annotated genes (at least compared to chromosomes), but all the "primary" ones should be in any GTF. I'm not even sure to what extent non-primary contigs (extra haplotypes, patches, etc.) are annotated and if there would be a GTF available to cover them. So I think you'll be fine with a GTF e.g. from Ensembl. If you are using VEP annotation of your variants, it might make sense to use the same Ensembl release for both the VEP cache and the GTF, so you get fully compatible gene & transcript IDs. We had this discussion a year ago, you might find it useful.

I know that @Louis Le Nézet has some scripts to remove the chr_ prefix in vcf for #C05G7UE94CT if I'm not mostaken

Thanks a lot for the quick reply @Thomas Danhorn! Your answer is very useful and well explained! Thanks for the suggestion @Maxime Garcia! I will definitely check the scripts! One question: In this case would be also possible to add the "chr" prefix to the BAM files instead of modifing all the GATK reference files, Which one is the preferable solution between modifying the BAM or all the GATK reference files?

In my opinion all references need to match

Thanks a lot for your help!

I don't think it matters all that much, it's more a question of the effort required. If you have dozens of BAMs, it might be easier to modify a couple of reference files, but if you just have a few BAMs (use reheader to modify), keeping the reference files intact has the advantage that you have fewer copies of reference files to keep track of ...

Thanks for all the suggestions!

Hi, we do effectively have some script to update bam and vcf file to a common reference and style (with or without chr prefix). For the moment we use it as a separate workflow (name

chrcheck

) before phaseimpute but we might in the future separate it from phaseimpute as a standalone small pipeline / subworkflow that could be used by other workflow.

Seems very interesting @Louis Le Nézet, thanks! I look forward for the sub-workflow. If it can be useful, I have a sub-workflow to make BAM files compatible with GATK requirements (correcting 'chr' string and @RG line)