Heads up:
We’re moving the GATK website, docs and forum to a new platform. Read the full story and breakdown of key changes on this blog.
Notice:
If you happen to see a question you know the answer to, please do chime in and help your fellow community members. We encourage our fourm members to be more involved, jump in and help out your fellow researchers with their questions. GATK forum is a community forum and helping each other with using GATK tools and research is the cornerstone of our success as a genomics research community.We appreciate your help!

Test-drive the GATK tools and Best Practices pipelines on Terra


Check out this blog post to learn how you can get started with GATK and try out the pipelines in preconfigured workspaces (with a user-friendly interface!) without having to install anything.

Calling variants in RNAseq

2»

Comments

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi @npriedig, as shown in the filtering summary the main reasons you are losing reads are MAPQ and duplicates, so it's a quality control problem. You mention a high mapping rate, but you should check what portion of your mapped reads are actually uniquely mapped, since those are the only reads that will be usable for variant discovery. Regarding duplicates, it's either a problem of true duplication, or the prep method you use generates the appearance of duplication (as in amplicon-based sequencing). If the latter, you can disable the DuplicateRead filer (using -drf) but you should evaluate results carefully if you do.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi @npriedig, sorry for the very late response. I thought I had answered this and closed the issue ticket prematurely.

    In a nutshell, you have two "problems" leading to filtering out of many reads: one is the high rate of duplication (41.90% of total as indicated in the filtering summary) and the other is the also high rate of reads with a MAPQ below 20 (35.78% of total) which HaplotypeCaller disregards as useless.

    Regarding the duplication, one thing to keep in mind is that for RNAseq the calculation of duplication rates is confounded by the possibility that some apparent duplicates are actually the result of higher levels of transcription. We still prefer to filter out duplicates, but some people choose to disable the duplicate filter (using -drf DuplicateRead).

    I would say the other problem, of low MAPQs, is more a concern. It indicates that even though it looks like your data is mapped, there are lots of reads that are not well mapped. If you followed the pre-processing procedure that we recommend, you will have converted all MAPQs of 255 (which are the unique mappers that we can use for variant discovery) to a flat 60. Those reads will automatically pass HC's filter. So I would advise looking at the distribution of MAPQs in your files (both the original file mapped by STAR and the one obtained after pre-processing to check if the problem was there from the start or whether it arose on the way.

    It's also worth checking what is the actual depth (DP) at most of your variant sites. It's possible that you don't really have a problem, if you still have sufficient coverage despite all the filtering to make confident calls.

  • kchang3kchang3 Member

    I have followed the GATK best practice guideline for preparing the BAM file for mutation calling, and noticed that majority of the mutation fall in 3'UTR region, which isn't something i notice in DNAseq. Is this a common in mutation calling in RNAseq data?

    Top mutation types in my data
    3'UTR 42.79%
    Missense_Mutation 16.36%
    Intron 12.62%
    Silent 7.86%

    Issue · Github
    by Sheila

    Issue Number
    820
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    vdauwera
  • uki156uki156 Member
    edited April 2016

    Hi all, I'm new here, sorry if this has been asked before, but from what I can see it hasn't been so far.
    I tried running this workflow and I get an error at the GATK Indel Realigner tool. Basically what the error says is:

    Error caching SAM record xxxx, which is usually caused by malformed SAM/BAM files in which multiple identical copies of a read are present.

    Now, I could skip this step as Indel Realigner is optional, but I'm guessing that further GATK tools in the workflow will give the same error, as this seems to be a problem with my aligned BAM file coming out of the STAR aligner.
    Basically, from what I could gather so far, my BAM file has multiple identical copies of a read, which I am not really sure what the problem is about. I thought the error could be from the Mapping Qualities which aren't proper, so I tried changing the parameter --outSAMmapqUnique in STAR from 255 to 60 as advised for these type of tasks, but I still get the same error.
    Also this only happens when I try this analysis with paired-end reads, single-end reads come through ok.
    Any ideas?

  • SheilaSheila Broad InstituteMember, Broadie admin

    @uki156
    Hi,

    What kind of data are you working with and how was it sequenced? Did you follow the Best Practices pre-processing steps?

    Thanks,
    Sheila

  • npriedignpriedig Member
    edited April 2016

    Hello @Geraldine_VdAuwera, thank you for getting back regarding my issue! What you wrote makes perfect sense and is very clear. I will look very carefully at the MAPQ values to see if I can improve the results.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @kchang3 Sorry for the late reply. This is not something I'm aware of, but we don't deal with RNAseq much so I really can't speak with any authority on the matter. My only thought is that this does smell a bit like an artifact, so maybe consider checking whether the distribution of annotation values and scores of variants that are falling in 3'UTR regions is in line with the rest or not. If there's something wrong with those variants, that might come up as a distinct profile on one or more dimensions.

  • JPEJPE MaltaMember

    What is the visualization software used in this tutorial? The one showing the read alignment?

  • SheilaSheila Broad InstituteMember, Broadie admin

    @JPE
    Hi,

    We use IGV, which is developed at the Broad.

    -Sheila

  • everestial007everestial007 GreensboroMember ✭✭

    I am getting some error while running variant filtration.
    WARN 22:37:17,483 Interpreter - ![0,2]: 'QD < 2.0;' undefined variable QD

    I have been trying to find the source of the error but no luck. Never had this kind of issues before with other vcf files when running a standard (best practices). Also, the QD field is indeed present in the raw.vcf file.
    Any suggestions is appreciated.

    Thanks,

  • SheilaSheila Broad InstituteMember, Broadie admin

    @everestial007
    Hi,

    Have a look at my response in this thread.

    -Sheila

  • everestial007everestial007 GreensboroMember ✭✭

    Thanks @Sheila !
    It really helped clear my concerns.

  • Hi, I'm new here. We are planning to make SNP calling in RNA seq data. We don't have a reference genome, we want to compare SNP between different species with different ploidy levels.
    We need to build a transcriptome reference or can we work with a related species?

    Any suggestions is really appreciated.

    Thank you very much.

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭
    edited October 2016

    Hi @Carolina_O,

    Our tools rely on reference-aligned sequences and calling variants against the reference. So depending on what type of variants you are interested in, e.g. those that vary from the transcriptome or those that vary between species, I will tentatively say that either should work for you. I suppose calling variants against the transcriptome may be simpler as alignment records would be kept minimal.

    Your project sounds complicated and we have some forum threads that should be of interest to you, e.g. here, here and here. Also, see this FAQ. I'm sure there are more discussions that would interest you and you should find them by searching the forum.

    Post edited by shlee on
  • Thank you very much Shlee!

  • vyellapavyellapa Member

    Curious if the haplotype caller canbe replaced with Mutect to identify low VAF variants?

  • SheilaSheila Broad InstituteMember, Broadie admin

    @vyellapa
    Hi,

    Have a look at this thread which should help.

    -Sheila

  • jealesjeales Member

    Hi,

    The STAR documentation strongly suggests supplying a GTF file of transcript annotations when building the genome index.
    https://github.com/alexdobin/STAR/blob/master/doc/STARmanual.pdf

    Is there any reason that step 1 does not supply a human GTF?
    I was hoping that if I supplied a GTF then the splice junctions used in the second pass would be more stable across multiple samples

    thanks
    James

    Issue · Github
    by Sheila

    Issue Number
    1365
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    vdauwera
  • JPEJPE MaltaMember

    Is it possible to annotate the resulting VCF with gene name and RS number (for human samples)? Which part of the analysis needs tweaking for this, the STAR alignment? Can you please provide an example?

  • ypnngaaypnngaa DCMember

    Hello,
    If I am also interested in phasing of those SNPs that called from this pipeline, can I just enable --emitRefConfidence GVCF? Or I just use the exactly same pipeline, use ReadsBackedPhasing?
    Thanks!

  • SheilaSheila Broad InstituteMember, Broadie admin

    @JPE
    Hi,

    You can use VariantAnnotator with --dbsnp, --resource, and --snpEffFile after calling variants with HaplotypeCaller.

    -Sheila

  • SheilaSheila Broad InstituteMember, Broadie admin

    @ypnngaa
    Hi,

    Yes, the latest version of HaplotypeCaller does phasing automatically, either in normal mode or in GVCF mode :smile:

    You only need to run ReadBackedPhasing if you want to merge MNPs.

    -Sheila

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @jeales I think this recommendation to use a GTF file came up after we wrote up our docs on this. Ultimately we defer to Alex Dobin for the most up to date way to run STAR since he is the author of that program. We'll try to make that more clear in the docs.

  • ypnngaaypnngaa DCMember
    edited October 2016

    Hi @Sheila ,
    Thanks for replying. However, in my case, it did not automatically perform phasing. I used this command to run on my RNA Seq data:
    java -jar GenomeAnalysisTK.jar -T HaplotypeCaller -R $genomeFasta -I $Input_BAM -dontUseSoftClippedBases -stand_call_conf 0 -stand_emit_conf 0 -o $Output_VCF
    In the output VCF, I did not see SNPs being phased. And I saw this message from GATK:
    HaplotypeCaller - Disabling physical phasing, which is supported only for reference-model confidence output
    Do you know why?
    If I want it to do phasing, what should I use?(I know it is RNAseq based, phasing can be conducted just locally, but I need that.) Now I just used seperated steps: first use UnifiedGenotype, then ReadBackedPhase. Do you have any suggestion?

    Thanks!

  • SheilaSheila Broad InstituteMember, Broadie admin

    @ypnngaa
    Hi,

    Oh, I see. I'm sorry I mislead you. It seems phasing is only supported in GVCF mode. Right now, we recommend HaplotypeCaller in normal mode only for RNA-seq, because we have not done testing for the GVCF workflow. However, some users have reported using the GVCF workflow for RNA-seq data.

    You can either try the GVCF workflow and thoroughly test your results, or run HaplotypeCaller then ReadBackedPhasing afterwards.

    -Sheila

  • J_RJ_R NCMember

    Hi,

    Yet another question, sorry. Is there more documentation unpacking why there's an RNA-specific recommendation of hard-filtering clustered SNPs?

    Thanks!

    Joanna

  • SheilaSheila Broad InstituteMember, Broadie admin

    @J_R
    Hi Joanna,

    I think this article is all we have for now. There may be more to come in the future.

    -Sheila

  • mjtivmjtiv Newark, DEMember
    edited January 2017

    Question about GATK 3.7 and RNA-Seq pipeline, is it safe to say the command "-stand_emit_conf 20.0" for STEP 6 has been removed from the pipeline commands? I got a warning this command was deprecated and so I removed it, otherwise the pipeline won't run.

    Just read your updates on the new version---- Which state you removed it.... My fault

    Post edited by mjtiv on
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin
    We'll update the article too, thanks for pointing this out.
  • everestial007everestial007 GreensboroMember ✭✭
    edited January 2017

    @Geraldine_VdAuwera @Sheila @shlee

    Regarding the filtering of VCF from RNAseq using SNPcluster:
    java -jar GenomeAnalysisTK.jar -T VariantFiltration -R hg_19.fasta -V input.vcf -window 35 -cluster 3 -filterName FS -filter "FS > 30.0" -filterName QD -filter "QD < 2.0" -o output.vcf

    Is there a method or functionality to prepare the bed file based on this SNP cluster?

    I would rather want to filter the reads from the bam file and run Denovo assembly and blast search to find other potential mapping areas for the reads on the genome.

    Since, this problem addresses paralogous alignment in which one set of reads (mostly containing reference allele) will be true alignment while the other would be paralogs; I was thinking if reads can be filtered in following manner.

    • select the site that have SNP clusters and prepare bed file.
    • A. filter reads based on that bed file.
    • B. Or, only select the reads that contains dominantly non-reference SNP cluster, suggesting that is the paralogous read. This way only paralogous alignment will be filtered.

    Issue · Github
    by Sheila

    Issue Number
    1713
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    sooheelee
  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    Hi @everestial007,

    If by paralogous alignments you are referring to reads that can multiply map within a given genome, then I think low or zero mapping scores (MAPQ) in a region would be the basis for your analysis workflow. I have to say that it's not clear to me how filtering on FisherStrand helps towards identifying regions with low/zero MAPQ reads. I do see that QualByDepth, being normalized by depth supporting a variant, could help but this metric also factors for quality. It's not clear then how quality would correlate with regions with reads with low or zero MAPQ scores. I'm just learning about these annotation metrics so perhaps I've misunderstood your approach.

    Can you please define what you mean by SNP cluster? We have some filters that target reads that are likely contaminant (GATK engine read filters here) as well as sites with clustered events (in MuTect2).

    Let me answer the parts of your questions that I can.

    Within Picard and GATK tools, we typically use interval lists more so than bed files. Here are some tools that can help you prepare an intervals list:

    To convert an intervals list to a bed file, please see the awk command in Tutorial#7156, right before the start of section 3.

    You might also be interested in Picard's MakeSitesOnlyVcf, to simplfy the starting VCF. Also, it appears Picard's FilterVcf may also do your FS/QD filtering.

    All this being said, to subset out reads to filter reads based on sites within a VCF you really don't need to create a bed/intervals file. You can simply supply the vcf with the -L parameter, e.g. -L file.vcf, in conjunction with a GATK tool, e.g. PrintReads. The -L option also accepts interval lists and bed files. You might want to double-check if all that is needed is overlap to subset reads or if you need to supply padding (--interval_padding) to cover the entirety of reads and perhaps change the merging approach (--interval_set_rule). These are defined in the GATK engine parameters that you can use with all GATK tools.

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    P.S. @everestial007. I've just confirmed with our engineer that you can print out all reads that touch the interval/site. No need to expand intervals.

  • mjtivmjtiv Newark, DEMember

    In your RNA-Seq pipeline it is mentioned known issues may occur with HaplotypeCaller with different alleles that are far from the 0.5 ratio which will influence ASE studies. The suggestion is to use a "DNA-aware" mode of the caller to fix this issue. Could you specifically explain more about this "DNA-aware" mode and how to perform it in HaplotypeCaller or if necessary the entire pipeline? Thank You

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @mjtiv What you're referring to is a hypothetical feature we were considering but never had the opportunity to build. The idea was to use a DNA normal to control for estimating baseline expectations.

  • JudyJudy ChinaMember

    I have the same problem with "MESSAGE: Error caching SAM record xxx, which is usually caused by malformed SAM/BAM files in which multiple identical copies of a read are present." during IndelRealigner, but I have tried adding command line -rf NotPrimaryAlignment , it doesn't work. Is there anybody having some solutions?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin
    @Judy please do NOT post the same comment multiple times. It clutters the forum and wastes our time.

    Have you validated your bam file?
  • kulbirsinghsandhukulbirsinghsandhu winnipegMember

    Hi everyone, I am doing trying to find variants using my RNAseq data (it was generated for expression analysis). I had 8 samples which were indexed and sequenced in a single lane. I received 16 FASTQ paired end files from 8 samples after sequencing. These 8 samples are derived from 2 vareities and 2 conditions. I am interested in finding variants among these 2 plant varaities. The fastq files were trimmed for adaptor sequences and only paired sequences were retained by using Trimmomatic. My question is that during the first and second alignment step with STAR, can i use all the fastq files for alignment at the same time? This will get me one output sam file to take into next steps upto base recalibration. Or should I align each read pair separately and then go upto base recalibration steps separetly for each read pair, and merge before variant calling step? Please suggest. Thanks, Kulbir

  • SheilaSheila Broad InstituteMember, Broadie admin

    @kulbirsinghsandhu
    Hi Kulbir,

    I think this thread will help.

    -Sheila

  • YogeshYogesh south koreaMember

    I am Trying to find SNP in the transcriptome. I have finished the Split'N'Trim and reassign mapping qualities step. after this can I do indel realingment using these command: But sorted bam file is needed, can I sort the bam file using samtool:

    /home/samtools/samtool sort dedup.bam >sorted.dedub.bam

    Generate intervals of interest from sample alignments:
    java -jar GenomeAnalysisTK.jar \
    -T RealignerTargetCreator \
    -nt 4 \
    -R refgenome.fa \
    -I **sorted.dedup.bam **\
    -o realign.intervals

    Realign (multiple sequence alignment):
    java -jar GenomeAnalysisTK.jar \
    -T IndelRealigner \
    -R refgenome.fa\
    -targetIntervals realign.intervals \
    -I sample1.sorted.dedup.bam \
    -o sample1.sorted.dedup.realigned.bam

    Fix mate pair info in BAM (PICARD):
    java -jar $PICARDDIR/picard.jar FixMateInformation \
    INPUT=sample1.sorted.dedup.realigned.bam \
    OUTPUT=sample1.sorted.dedup.realigned.fixmate.bam \
    SO=coordinate \
    CREATE_INDEX=true

  • GovardhanGovardhan IndiaMember

    Hi,

    I performed variant calling using RNA-seq in our cohort as we had poor quality exome seq data and thanks to GATK for this pipeline I identified commonly identified hot-spot mutations.

    We are planning to include these findings in upcoming publication and I have question about cutting this pipeline.

    Should I just cite three GATK publication??

    Also, did anyone published anything using same pipeline??

    Thanks in advance.

    Cheers,
    Gova

  • SheilaSheila Broad InstituteMember, Broadie admin

    @Govardhan
    Hi Gova,

    We recommend citing both the Best Practices pipeline (from the GATK site) and one of the GATK papers.

    -Sheila

  • mjtivmjtiv Newark, DEMember

    With variant calling step (HaplotypeCaller), how is strand direction taken into account? The reason I am asking is because I am trying to compare RNA-seq variant data produced using this pipeline to genotyping data that is always (+) strand.

  • everestial007everestial007 GreensboroMember ✭✭

    @mjtiv

    I think the variants are always reported in context of the reference genome you supplied. The RNAseq transcripts may come from either + or - strand but, variants are always reported in context of the ref allele.

  • everestial007everestial007 GreensboroMember ✭✭

    @mjtiv
    Hope you already know this now. But, just to verify that assertion I made earlier: https://groups.google.com/forum/#!topic/rna-star/ehHFgV_2cjU
    See the response on similar problem, from Alex Dobin the author of STAR-Rna.

  • Hi, has this pipeline been tested with the HISAT2 aligner? If so, how did it compare to the results with STAR?

  • SheilaSheila Broad InstituteMember, Broadie admin

    @lshepard
    Hi,

    According to this article, "For RNA-seq, we evaluated all the major software packages that are specialized in RNAseq alignment, and we found that we were able to achieve the highest sensitivity to both SNPs and, importantly, indels, using STAR aligner. ". You can find more details under "1. Mapping to the reference"

    -Sheila

  • lshepardlshepard Member
    edited August 2017

    @Sheila said:
    @lshepard
    Hi,

    According to this article, "For RNA-seq, we evaluated all the major software packages that are specialized in RNAseq alignment, and we found that we were able to achieve the highest sensitivity to both SNPs and, importantly, indels, using STAR aligner. ". You can find more details under "1. Mapping to the reference"

    -Sheila

    @Sheila
    Yes, I noticed that, but I am it is not very clear to me when this note was added to the page as it was first created in 2014 when HISAT2 was not available yet. It says it was last updated in 2017, but were the newer aligners tested as well?

    Thanks for the clarification!

    Issue · Github
    by Sheila

    Issue Number
    2461
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    sooheelee
  • SheilaSheila Broad InstituteMember, Broadie admin

    @lshepard
    Hi,

    I am not sure if there has been any testing that recently. I will check with the team and get back to you.

    Of course, STAR is what we recommend, but you are free to use whatever aligner you have determined is best for your needs :smile:

    -Sheila

  • @Sheila

    Thank you! I appreciate you check-in with the team!

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    Hi @lshepard,

    I believe testing was done between TopHat and Star and the decision to use Star had in part to do with its output of compatible metadata, e.g. CIGAR string representations.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    That's right, a comparison was made on several packages including TopHat and Star; criteria included sensitivity metrics and interoperability of results (well-behaved output format). Wrt sensitivity we found that Star gave the best overal results across SNPs and indels, though iirc other packages may have provided slightly better results over SNPs alone. That was several years ago though, and we have not reevaluated this since as RNAseq has not been a development priority. This is likely to change in the near future however.

  • lshepardlshepard Member
    edited September 2017

    @Geraldine_VdAuwera
    Sounds good, thank you for giving the details of the analysis performed! I might do some of these comparisons myself in the near future.

  • krdavkrdav Member, Broadie

    @lshepard
    Somebody else has already done a similar comparisons that indeed does claim better results using TopHat2:
    https://wellcomeopenresearch.org/articles/2-6/v2

  • Hi Geraldine, I am planning to do variant calling in RNA-seq data but haven't prepared the sample library yet. Could you please give me some advice on the minimum sequencing depth required to find meaningful true variants?

  • lshepardlshepard Member
    edited September 2017

    @krdav said:
    @lshepard
    Somebody else has already done a similar comparisons that indeed does claim better results using TopHat2:
    https://wellcomeopenresearch.org/articles/2-6/v2

    @krdav
    Interesting, thanks for pointing it out. I will need to check the details, but if that is indeed the case then I would expect potentially better results with HISAT2. Thanks again!

  • SheilaSheila Broad InstituteMember, Broadie admin

    @uff
    Hi,

    We don't have any specific recommendations, but GATK tools are tested on 30X data.

    -Sheila

  • Hi @Geraldine_VdAuwera and @Sheila,

    It has been several years since this article was originally was posted. You described this possibility: "We also plan to add functionality to process DNAseq and RNAseq data from the same samples simultaneously"

    Is HaplotypeCaller yet capable of using WGS and RNA-seq data simultaneously to improve SNP/Indel calling?

    Thanks!
    Dan

    Issue · Github
    by Sheila

    Issue Number
    2511
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    chandrans
  • SheilaSheila Broad InstituteMember, Broadie admin

    @dsap
    Hi Dan,

    I don't think there has been much done in that area, so it is not a good idea to run on both DNA and RNA together. Let me confirm with the team and get back to you.

    -Sheila

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Unfortunately we haven't implemented that functionality, sorry.

  • jgouldjgould GouldMember ✭✭

    Is there a WDL available that implements this workflow?

    Issue · Github
    by Sheila

    Issue Number
    2608
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    vdauwera
  • SheilaSheila Broad InstituteMember, Broadie admin

    @jgould
    Hi,

    I am checking with the team and will get back to you.

    -Sheila

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @jgould We currently don't have a WDL for variant calling on RNAseq; it's on the roadmap of workflows we aim to provide and support in early 2018.

  • lidialidia Member

    Hello GATK team I am new in the use of this toolkit and therefore I have some questions. I am working with RNA paired end read sequences. My goal is to find SNPs in 2 strains of trout to be able to compare them, from each strain I have 8 fish sampled, each fish with 6 different tissues. I already went through the best practices workflow and I have several questions about the STAR steps. My first steps are the following:

    picard CreateSequenceDictionary R= genome.fna O= genome.dict

    samtools faidx ReferenceGenome/genome.fna

    STAR --runMode genomeGenerate --genomeDir GenomeDir –genomeFastaFiles ReferenceGenome/genome.fna --sjdbGTFfile ReferenceGenome/Annotations.gff --limitGenomeGenerateRAM 140000000000 --runThreadN 8

    STAR --runMode alignReads --genomeDir GenomeDir --readFilesCommand zcat --readFilesIn Forelle/${R1} Forelle/${R2} --outFileNamePrefix ${sample}_ --runThreadN 8 /

    I understand that I have to do the previous step for every tissue independently. Therefore I would be generating 96 SJ.out.tab files.

    For the second pass how can I create the genome indexes if I have 96 files of spliced junctions? or do I have to repeat this process independently 96 times so that it will be like creating 96 genome directories?

    STAR --runMode genomeGenerate --genomeDir GenomeDir2 --genomeFastaFiles ReferenceGenome/genome.fna --sjdbFileChrStartEnd ${sample}_SJ.out.tab --sjdbOverhang 64 --limitGenomeGenerateRAM 140000000000 --runThreadN 8 /

    STAR --genomeDir GenomeDir2 --readFilesCommand zcat --readFilesIn Forelle/$R1 Forelle/$R2 --outFileNamePrefix ${sample}2pass --runThreadN 8 /

    Thanks!

  • SheilaSheila Broad InstituteMember, Broadie admin

    @lidia
    Hi,

    I am not sure I understand what you mean by "For the second pass how can I create the genome indexes if I have 96 files of spliced junctions? or do I have to repeat this process independently 96 times so that it will be like creating 96 genome directories?" What do you mean by "genome indexes"? Are you aligning to one reference genome?

    Thanks,
    Sheila

  • lidialidia Member
    edited January 2018

    Hi Sheila,

    Yes, I am aligning to one reference genome. Therefore, since I understand that the alignment has to be done tissue by tissue, at the end of the first pass I have 96 SJ.out.tab files. But I am unsure about how to include them for the second pass.

  • SheilaSheila Broad InstituteMember, Broadie admin

    @lidia
    Hi,

    Ah, I see. Thanks for the clarification. We have an example of how to run the alignment step here. I think you will have to re-run for 96 files. You may want to ask the STAR team for any other questions about the usage, as we cannot help much with that.

    -Sheila

    P.S. For your samples, are you trying to determine differences in tissues? I am assuming so. You may also find this thread helpful.

  • Hi GATK Team,

    It seems like you guys didn't use GTF annotations in the STAR index mapping. Would you guys recommend using GTF for mapping RNA to reference using STAR? Some of our co-workers suggested that it might be an option. Would that improve the variant calling quality?

    Thanks

  • bhaasbhaas Broad InstituteMember, Broadie

    Is there an updated version of this protocol for use with gatk-4, as some of the options in gatk have changed.

  • SheilaSheila Broad InstituteMember, Broadie admin

    @flynnchen
    Hi,

    I think this thread may be helpful.

    -Sheila

  • SheilaSheila Broad InstituteMember, Broadie admin

    @bhaas
    Hi,

    The Best Practices are in development for GATK4, so we don't have full documentation for it yet.

    -Sheila

  • VrushaliVrushali Member

    Hi,

    The tools RealignerTargetCreator and IndelRealigner are obsolete in GATK4. Are there any replacements for these tools in GATK4?

    Also, there are two option -BQSR (used with PrintReads) and -stand_emit_conf (used with HaplotypeCaller) that are not available with GATK4. What should be done in these cases?

  • SheilaSheila Broad InstituteMember, Broadie admin

    @Vrushali
    Hi,

    Have a look at this thread. However, there are plans to port those tool now, and you can keep track of it here.

    For BQSR, have a look at this thread. For stand_emit_conf, have a look at this blog post to find out why it is gone.

    -Sheila

  • mjtivmjtiv Newark, DEMember
    edited February 2018

    Question about pipeline, I just got into a conversation with a lab mate about using a gtf file in an RNA analysis pipelines and realized that GATK's Best Practices for RNA-Seq does NOT appear to use a reference gtf file in the analysis (slight heart attack thinking need to re-run everything). Could you expand on this gtf free option? Also, I did look back at my notes and see I did use gtf file to help create the STAR genome index reference file (hopefully doesn't mess things ups).

    Any feedback on this topic would be greatly appreciated and fingers crossed its not bad news.

  • SheilaSheila Broad InstituteMember, Broadie admin

    @mjtiv
    Hi,

    Have a look at this thread.

    -Sheila

  • lfalllfall Member

    Hello,
    I am trying to call SNPs in my RNA-seq data using the GATK pipeline, but I am running into the same error over and over. The alignments were done using STAR and I was able to make it through the workflow steps 1 (mapping to reference (STAR)) and two (Add read groups, sort, mark duplicates, create index (Picard)). When I get to the Split'N'Trim step, I get an error. I have tried changing and removing just about everything I can think of in the command but I still get this error saying that the argument is wrong. I've checked so many times on the GATK best practices to see if my scripts are wrong, but they seem to match what is written on the workflow. I'm very sorry if this question has an obvious answer, I have tried so many things and I can't seem to find the answer anywhere! Any help is very much appreciated :) Thank you!

    This is the script I am running:

    !/bin/bash

    SBATCH -N 1

    SBATCH --ntasks 2

    SBATCH --mem 24G

    SBATCH --mail-type ALL

    SBATCH -J split.n.trim

    module load GATK/3.8-0-Java-1.8.0_121
    module load Java/1.8.0_152
    module load R/3.2.5-IGB-gcc-4.9.4
    module load Python/2.7.13-IGB-gcc-4.9.4

    java -jar GenomeAnalysisTK.jar \
    -T SplitNCigarReads \
    -R /home/a-m/lfall/Fusarium_rnaseq/data/genome/NCBI/GCF_000240135.3_ASM24013v3_genomic.fna \
    -I BHS.dedupped.bam \
    -o BHS.split.bam \
    -rf ReassignOneMappingQuality -RMQF 255 -RMQT 60 -U ALLOW_N_CIGAR_READS

    This is the error I get after running this script:

    [[email protected] rg.duplicates.index]$ more slurm-622962.out
    To execute GATK run: java -jar $EBROOTGATK/GenomeAnalysisTK.jar

    The following have been reloaded with a version change:
    1) Java/1.8.0_121 => Java/1.8.0_152

    The following have been reloaded with a version change:
    1) Java/1.8.0_152 => Java/1.8.0_121

    INFO 16:00:00,666 HelpFormatter - ----------------------------------------------------------------------------------
    INFO 16:00:00,685 HelpFormatter - The Genome Analysis Toolkit (GATK) v3.8-0-ge9d806836, Compiled 2017/07/28 21:26:50
    INFO 16:00:00,685 HelpFormatter - Copyright (c) 2010-2016 The Broad Institute
    INFO 16:00:00,685 HelpFormatter - For support and documentation go to https://software.broadinstitute.org/gatk
    INFO 16:00:00,685 HelpFormatter - [Thu May 10 16:00:00 CDT 2018] Executing on Linux 3.10.0-693.11.1.el7.x86_64 amd64
    INFO 16:00:00,685 HelpFormatter - Java HotSpot(TM) 64-Bit Server VM 1.8.0_121-b13
    INFO 16:00:00,688 HelpFormatter - Program Args: -T SplitNCigarReads -R /home/a-m/lfall/Fusarium_rnaseq/data/genome/NCBI/GCF_000240135.3_ASM24013v3_geno
    mic.fna -I BHS.dedupped.bam -o BHS.split.bam -rf ReassignOneMappingQuality -RMQF 255 -RMQT 60 -U ALLOW_N_CIGAR_READS
    INFO 16:00:00,694 HelpFormatter - Executing as [email protected] on Linux 3.10.0-693.11.1.el7.x86_64 amd64; Java HotSpot(TM) 64-Bit Server VM 1.8.0_12
    1-b13.
    INFO 16:00:00,699 HelpFormatter - Date/Time: 2018/05/10 16:00:00
    INFO 16:00:00,699 HelpFormatter - ----------------------------------------------------------------------------------
    INFO 16:00:00,699 HelpFormatter - ----------------------------------------------------------------------------------
    ERROR StatusLogger Unable to create class org.apache.logging.log4j.core.impl.Log4jContextFactory specified in jar:file:/home/a-m/lfall/Fusarium_rnaseq/r
    esults/NCBI_GATK/start.from.sam/rg.duplicates.index/GenomeAnalysisTK.jar!/META-INF/log4j-provider.properties
    ERROR StatusLogger Log4j2 could not find a logging implementation. Please add log4j-core to the classpath. Using SimpleLogger to log to the console...
    INFO 16:00:00,997 GenomeAnalysisEngine - Deflater: IntelDeflater
    INFO 16:00:00,997 GenomeAnalysisEngine - Inflater: IntelInflater
    INFO 16:00:00,997 GenomeAnalysisEngine - Strictness is SILENT

    ERROR ------------------------------------------------------------------------------------------
    ERROR A USER ERROR has occurred (version 3.8-0-ge9d806836):
    ERROR
    ERROR This means that one or more arguments or inputs in your command are incorrect.
    ERROR The error message below tells you what is the problem.
    ERROR
    ERROR If the problem is an invalid argument, please check the online documentation guide
    ERROR (or rerun your command with --help) to view allowable command-line arguments for this tool.
    ERROR
    ERROR Visit our website and forum for extensive documentation and answers to
    ERROR commonly asked questions https://software.broadinstitute.org/gatk
    ERROR
    ERROR Please do NOT post this error to the GATK forum unless you have really tried to fix it yourself.
    ERROR
    ERROR MESSAGE: Invalid command line: Failed to load reference dictionary
    ERROR ------------------------------------------------------------------------------------------
  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    Hi @lfall,

    The message:

    ERROR MESSAGE: Invalid command line: Failed to load reference dictionary

    Indicates the tool is unable to access the dictionary associated with the reference. Your reference is

    -R /home/a-m/lfall/Fusarium_rnaseq/data/genome/NCBI/GCF_000240135.3_ASM24013v3_genomic.fna

    So here are the things I would try:

    1. Generate a reference dictionary with Picard CreateSequenceDictionary in the same directory as the reference.
    2. It's possible the tool is having an issue with the .fna extension of your reference. Typical reference FASTA extensions are .fasta or .fa. Try switching this out, assuming your reference is actually in simple FASTA format. Otherwise, try to convert your FNA file to FASTA.
  • after creating an index with only 1.sj.out file. I realized that my sam files are missing columns. How can I create the index with multiple sj.out.tab files ?

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    Hi @kibrahim90, it's not clear what you are asking. Perhaps you will find SAM specifications here helpful.

  • kibrahim90kibrahim90 Member
    edited May 2018

    @shlee said:
    Hi @kibrahim90, it's not clear what you are asking. Perhaps you will find SAM specifications here helpful.

    how do I run star (aligner) to create an index with multiple sj.out.tab files is what im asking. Its not in the manual. I checked.

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    @kibrahim90, Star is not a GATK/Picard tool and therefore we do not support it. Please try posting your question to Biostars or SeqAnswers.

  • CarolineJCarolineJ Member
    edited July 2018

    Hi team -

    Do you have any suggestions for the variant calling step (#6) and the hard filtering step (#7) for multisample RNAseq experiment aligned to a de novo transcriptome? Also, is there a still a benefit to using HaplotypeCaller over UnifiedGenotyper for cases where the reference is a de novo assembled transcriptome?

    1. GATK best practices for variant calling for a single sample/align to ref. genome
    java -jar GenomeAnalysisTK.jar -T HaplotypeCaller -R ref.fasta -I input.bam -dontUseSoftClippedBases -stand_call_conf 20.0 -o output.vcf
    

    My "best guess" at customizing GATK best practices for a multisample experiment with alignment to a de novo transcriptome

    java -jar GenomeAnalysisTK.jar -T HaplotypeCaller -R de.novo.transcriptome.fasta -I input_sample1.bam input_sample2 -stand_call_conf 20.0 -o output.vcf
    
    1. GATK best practices for hard filtering
    java -jar GenomeAnalysisTK.jar -T VariantFiltration -R de.novo.transcriptome.fasta -V input.vcf -window 35 -cluster 3 -filterName FS -filter "FS > 30.0" -filterName QD -filter "QD < 2.0" -o output.vcf 
    

    My best guess is that the hard filtering settings would be identical for a multisample/denovo transcriptome ref design.

  • SheilaSheila Broad InstituteMember, Broadie admin

    @CarolineJ
    Hi,

    Unfortunately, our team has not dedicated much time to the RNA pipeline, but we do have some recommendations in the github repo here. I think the team is now working on better Best Practices for RNA-seq, so keep an eye out for those.

    -Sheila

  • Hello,

    When I run the command

    mkdir runDir
    cd runDir
    STAR --genomeDir /scratch/snyder/y/yin168/gatk/test/genomeDir --readFilesIn /scratch/snyder/y/yin168/gatk/test/400_R1.fastq /scratch/snyder/y/yin168/gatk/test/400_R2.fastq --runThreadN 5

    I always got an error:

    /var/spool/torque/mom_priv/jobs/1411989.snyder-adm.rcac.purdue.edu.SC: line 18: 34161 Segmentation fault

    (sometimes it is line 16, instead of line 18)

    How should I resolve this issue?

    Thanks a lot!

  • Hello,

    In the gatk RNAseq variant calling pipeline, there is a step for the 2-pass STAR as below:

    3) For the 2-pass STAR, a new index is then created using splice junction information contained in the file SJ.out.tab from the first pass:

    genomeDir=/path/to/hg19_2pass mkdir $genomeDir STAR --runMode genomeGenerate --genomeDir $genomeDir --genomeFastaFiles hg19.fa --sjdbFileChrStartEnd /path/to/1pass/SJ.out.tab --sjdbOverhang 75 --runThreadN

    Here, --sjdbOverhang is set to 75. However, the default value of --sjdbOverhang is 100. Why is it set to 75 instead of 100? Does that make a difference? What value should I set for my analysis.

    Thanks a lot!

  • XiaoshenYinXiaoshenYin Member
    edited November 2018

    Hello,

    When running the gatk RNAseq variant calling pipeline, I am not quite sure about whether I should spend time doing 4. Indel Realignment and 5. Base Recalibration.

    For 4 - How significant the effect would be if I ignore this step? It says "we only recommend performing the realignment step if you have compute and time to spare (or if it’s important not to miss any potential indels)". I am calling SNPs from RNAseq, so how should I evaluate/determine whether it is important not to miss any potential indels? Could you list potential cases in which it is important not to miss any potential indels (or in which it is okay to miss potential indels)?

    Also, I see a post saying that, if a haplotype assembly is done (for example, run HaplotypeCaller), then it is not necessary to run Indel Realignment. Is this true? Since I am going to do 6. Variant calling by running HaplotypeCaller, does this mean I do not need to run Indel Realignment at all?

    For 5 - How big/marginal the effect would be if I just skip 5. Base Recalibration? When would the qualities have systematic error modes, as mentioned in the explanation for this step? Also, since --known-sites is a parameter I would need to provide in Base Recalibration, but I do not have a "known sites" resource file. Does it still make sense to do this step? If so, how should I generate/prepare a "known sites" resource file?

    Thanks.

    Post edited by XiaoshenYin on
  • Hello,

    I use the command to filter called variants according to the pipeline:

    java -jar ~/bin/GenomeAnalysisTK-3.8-0-ge9d806836/GenomeAnalysisTK.jar -T VariantFiltration -R ../sl_genome/gPmar100.fa.txt -V 400.vcf -window 35 -cluster 3 -filterName FS -filter "FS > 30.0" -filterName QD -filter "QD < 2.0" -o 400_filter_5.vcf --setFilteredGtToNocall

    In the filter column, I see "PASS", "snpcluster" and "QD". I would assume "QD" indicates that a specific locus has a poor call (i.e. with QD<2.0) and thus this variant should be filtered. Even though I add "--setFilteredGtToNocall", I still see that the genotype is "0/1". However, the genotype of this variant should be "./." because I use "--setFilteredGtToNocall". I try different options even with the latest gatk 4, but it still does not work. How can I make --setFilteredGtToNocall work and filter out the poor variants or set their genotypes to ./.?

    Also, what is the difference between gatk-package-4.0.11.0-local.jar and gatk-package-4.0.11.0-spark.jar? The gatk-package-4.0.11.0-spark.jar does not work on my cluster, so is it okay if I use gatk-package-4.0.11.0-local.jar?

    Thanks.

  • Hello,

    I am calling variants from RNAseq data using this gatk pipeline. I did what the pipeline suggest exactly and did not set ploidy for HaplotypeCaller. I got both "." and "./." in the vcf file. "./." stands for a diploid site at which the genotype cannot be called and "." stands for a haploid site at which the genotype cannot be called. Why do I get both missing haploid and diploid sites? To be more specific, why do I get ".", haploid site at which the genotype cannot be called? Should I keep or remove "." in analysis?

    Thanks.

  • TomWillDoTomWillDo Member
    Hi there,

    Sorry if this has been answered previously, but I can't seem to find it. Previous SNP/Indel hard filters (to serve in place of VQSR) include filters for SOR, but they seem to be missing from the RNA-seq pipeline. If someone could point me towards why they haven't been included, I would be appreciative!
  • hamaorhamaor Member
    Thanks a lot for this guide!
    I used this pipeline to analyze 6 samples of RNA-Seq (3 replicate of treatment and 3 replicate of control).
    after generating the filtered vcf, I'm using the GenotypeGVCFs to join treatment replicates together and control replicates together. (and I might want to join the whole dataset).
    Technically, 3 vcf files each 1.7G was joined into single 24MB file.
    is that reduction normal for the algorithm? can I make any validation to make sure I'm on the safe side?

    thanks
    Maor
  • encoded79encoded79 Member
    edited January 23
    Hello there

    I used your mutation analysis pipeline for RNA-seq data. However, when i compare the BAMs , pre and post CIGAR hard clipping on the iGV , i still see overlapping reads for intronic regions. Basically, no clipping seems to have been done. This is also reflected on the VCF files, where nearly half of variants are in fact intronic. I appreciate it if you could elaborate on the possible reasons for this. Thanks a lot
  • encoded79encoded79 Member
    edited January 28
  • XiaoshenYinXiaoshenYin Member

    Hello @shlee ,

    Step 3 in the pipeline for calling variants from RNAseq recommends:

    java -jar GenomeAnalysisTK.jar -T SplitNCigarReads -R ref.fasta -I dedupped.bam -o split.bam -rf ReassignOneMappingQuality -RMQF 255 -RMQT 60 -U ALLOW_N_CIGAR_READS

    However, when I try gatk4.0, I did not find "-rf ReassignOneMappingQuality -RMQF 255 -RMQT" or "-U ALLOW_N_CIGAR_READS" in optional arguments, although they are present in gatk3.8. Is it because updates on gatk4.0 make it unnecessary to set those two parameters? In this case, can I use gatk 4.0 without setting those two arguments to run the step of Split'N'Trim and reassign mapping qualities?

    Thanks.

  • seven32seven32 Member
    What is the correct syntax for calling HaplotypeCaller for RNA-seq data for GATK 4?
  • AlFatlawiAlFatlawi Member
    @Geraldine_VdAuwera
    Hi, regarding your comment in Oct 2014,
    "At this time the RNAseq pipeline is only meant for single-sample processing, yes. We have not finished testing how the tools perform on RNAseq for multiple samples."

    Please, is there any update in GATK 4?
    I have many samples (RNA Seq) and I would like to merge them in a VCF file.

    Thanks in advance.
    Ali
  • AlFatlawiAlFatlawi Member
    Hi,
    At 7. Variant filtering:
    I see you recommended the following:
    --filterExpression "FS > 30.0" -filterName QD --filterExpression "QD < 2.0"

    Please, are these for SNPs, Indel or both?
    and, what is your opinion if I use the same filterexpression on both SNPs and Indel?
    Regards
    Ali
  • AlFatlawiAlFatlawi Member
    any answer, please!
  • zjwangzjwang ChinaMember
    > @XiaoshenYin said:
    > Hello @shlee ,
    >
    > Step 3 in the pipeline for calling variants from RNAseq recommends:
    >
    > java -jar GenomeAnalysisTK.jar -T SplitNCigarReads -R ref.fasta -I dedupped.bam -o split.bam -rf ReassignOneMappingQuality -RMQF 255 -RMQT 60 -U ALLOW_N_CIGAR_READS
    >
    > However, when I try gatk4.0, I did not find "-rf ReassignOneMappingQuality -RMQF 255 -RMQT" or "-U ALLOW_N_CIGAR_READS" in optional arguments, although they are present in gatk3.8. Is it because updates on gatk4.0 make it unnecessary to set those two parameters? In this case, can I use gatk 4.0 without setting those two arguments to run the step of Split'N'Trim and reassign mapping qualities?
    >
    > Thanks.
    Hello, I have the same question as you, Have you solved this problem? Why do not find "-rf ReassignOneMappingQuality -RMQF 255 -RMQT" or "-U ALLOW_N_CIGAR_READS" in optional arguments
  • XiaoshenYinXiaoshenYin Member
    edited October 31

    Hello, how should I properly cite the pipeline of Calling variants in RNAseq in a manuscript?

  • XiaoshenYinXiaoshenYin Member
    edited November 5

    @zjwang Hello, I ended up using the older version, if I remember correctly.

Sign In or Register to comment.