RNA-Seq Variant Calling
I have RNA-seq data from two separate projects:
== Project 1==
- We received an RNA sample for a particular antibody strain. Over 150M, PE150 reads were generated.
== Project 2 ==
- We prepared another library from the same RNA sample. Over 24M, PE300 reads were generated.
I used the The GATK Best Practices for variant calling on RNAseq pipeline to analyze both projects. In Project 1, the VCF file reported no variants. In Project 2, the VCF file reported 4 variants.
What's an explanation for the difference in the VCF reports when the same GATK pipeline was used for two libraries prepared from the same sample?
The commands I used: https://pastebin.com/zh07Dvmr. Note: I ran the pipeline using Snakemake.
Update: I posted my question to the GATK forums.
• 983 views
•
link
written
by
user5212 •
0 answers
No answers yet.
Log in to answer this question.
More posts like this
-
Read group addition order in RNA-seq short variant discovery (SNPs + Indels)
written by iamsmor •
-
best pipeline/practice for analyzing WES for germline and somatic variant calling
written by epigene
-
Using GATK MarkDuplicates for targeted sequencing data
written by asalimih •
-
When to use .vcf or .gvcf files from GATK HaplotypeCaller?
written by vitor •
-
VCF files never get annotated with PG tag for CalculateGenotypePosteriors
written by lait
-
RNA-SEQ SNP calling
written by babji2318 •
-
How to combine a multi sample VCF from multiple sample VCF and interpret a multi sample vcf of bact…
written by S AR
-
Too many intronic variants in RNASeq, is this ok?
written by Sharon
-
Suggested pipelines for finding somatic mutations using RNAseq in normal and tumor cells
written by Sharon
-
Why the Allele frequency changed when using different variant calling tools?
written by murali
Were the sequencing runs from the same RNA extraction? Different extractions, same individual / cell line? Or is the same strain, but different individuals / cell lines?
Open the bam files on IGV and carefully examine the variants, you may spot something missed by the analysis pipeline, or at least find a plausible explanation. Variant calling with RNAseq is substantially less precise than with DNaseq, and issues like unequal coverage due to differential expression and differential allele expression may cause this.
I followed up with the lab techs who did the library prep. They said sequencing runs were from the same RNA extraction. The same RNA sample was used to prepare two different libraries. In the first project, a library was prepared and sequenced at PE150. In the second project, a different library was prepared and sequenced at PE300. I updated the original post.
Small nitpick. There are no PE150/PE300 libraries. Two different libraries were made (perhaps with different size inserts, may be worth asking). They were then sequenced at 150 bp and 300 bp.
If you are only seeing 0 or 4 variants per sample, something is wrong. There should be at thousands (at least for mammalian genomes). Are these filtered? How are you filtering them?
I mostly followed the GATK best practices. The raw reads were trimmed (>Q30) using Cutadapt and mapped to the transcript sequence using STAR (GATK best practices). I used Picard to mark duplicates.
I used GATK HaplotypeCaller with the
--dontUseSoftClippedBasesto hard clip any reads that could potentially call false variants due to overhanging regions. I applied a GATK filter to identify variants with a >Q20 and a Fisher Strand (FS) > 30 (Sample_filtered.vcf).You also ran
VariantFiltrationandFastaAlternateReferenceMaker. How big were the original VCFs fromHaplotypeCaller?In Project 1, in the original VCF still no variants were reported. In Project 2, in the original VCF the same four variants were reported. Maybe the
stand_call_confflag I use in theHaplotypeCallercommand is too stringent? Thestand_call_confflag uses the "minimum phred-scaled confidence threshold".No. 20 is not very stringent.
Can you share the command lines you used?
The commands I used: https://pastebin.com/zh07Dvmr.
What's the quality of your data? What do your fastqc, alignment, and picard tools QC metrics look like?
d-cameron is right, start by doing a basic sanity check. Is it possible all your reads are garbage, and nothing aligned? Is it possible you aligned to the wrong reference?