SEQanswers

Go Back   SEQanswers > Applications Forums > RNA Sequencing



Similar Threads
Thread Thread Starter Forum Replies Last Post
FastQC analysis pari_89 Bioinformatics 10 05-18-2013 02:55 AM
FASTQC help arkal Bioinformatics 1 10-17-2012 12:31 AM
Fastqc report melNGS Bioinformatics 3 07-24-2012 12:10 AM
Where is FastQC? sklages General 10 02-05-2012 11:46 PM
fastQC papori RNA Sequencing 3 02-04-2012 01:48 PM

Reply
 
Thread Tools
Old 08-05-2014, 05:59 AM   #1
Isa0984
Member
 
Location: Germany

Join Date: Aug 2014
Posts: 15
Default FastQC

Hello,
we've done a RNA-Seq analysis (Illumina HiSeq2000, 50 bp, paired end) and I have checked the quality with FastQC. There raised some questions for me:

1. Is FastQC as quality check alright for paired end reads?
2. The program gives for all of my four samples a fail for the per base quality of the second read (only the last three bases show lower quartile less than 5 or a median less than 20). Is there a logical explanation?
3. The sequencing and mapping was done by a company and they told us they trimmed the adapters. But I get a fail in the dublication level and if you look at the overrepresented sequences I can see only primer or adapter sequences. Did they a bad job?

Thanks for your help! Isabelle
Isa0984 is offline   Reply With Quote
Old 08-05-2014, 06:14 AM   #2
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,060
Default

FastQC is appropriate for QC of PE reads.

It would be better if you post screenshots/images of the FastQC results instead of just descriptions. Having something marked as "fail" does not automatically fail the entire sample. It is possible that the analysis done by your provider may not have removed all adapter dimers etc.
GenoMax is offline   Reply With Quote
Old 08-06-2014, 01:23 AM   #3
Isa0984
Member
 
Location: Germany

Join Date: Aug 2014
Posts: 15
Default

Ok, here are the images of FastQC...
Attached Images
File Type: png per_base_quality-1.png (8.6 KB, 14 views)
File Type: png per_base_quality-2.png (8.9 KB, 12 views)
Isa0984 is offline   Reply With Quote
Old 08-06-2014, 01:34 AM   #4
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Read 2 often has decreased quality at its 3' end. A bit of trimming can easily get rid of that.

BTW, they likely sent you untrimmed sequences and aligned trimmed sequences, which is why fastQC is telling you that the raw sequences still have adapter contamination.

Also, a fail on duplication level is pretty much expected for RNAseq data (that test is really only meant for whole-genome sequencing).
dpryan is offline   Reply With Quote
Old 08-06-2014, 03:15 AM   #5
Isa0984
Member
 
Location: Germany

Join Date: Aug 2014
Posts: 15
Default

Thanks a lot, that looks for me that the quality check makes not really sense then, its more or less good for the per base quality... ?
Isa0984 is offline   Reply With Quote
Old 08-06-2014, 03:28 AM   #6
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Yeah, just do a bit of quality/adapter trimming (e.g., with trimmomatic or trim_galore) and you should be good to go.
dpryan is offline   Reply With Quote
Old 08-06-2014, 03:36 AM   #7
Isa0984
Member
 
Location: Germany

Join Date: Aug 2014
Posts: 15
Default

But can I be shure that the company used trimmed data for mapping? Maybe they didnt, how can I check this?
Isa0984 is offline   Reply With Quote
Old 08-06-2014, 03:43 AM   #8
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Just look at the read lengths in the BAM file:

Code:
samtools view some_file.bam | cut -f 10 | awk '{print length($1)}' | uniq | sort | uniq
If they trimmed the reads prior to alignment, you should get more than one value.
dpryan is offline   Reply With Quote
Old 08-06-2014, 04:08 AM   #9
Isa0984
Member
 
Location: Germany

Join Date: Aug 2014
Posts: 15
Default

I will, but unfortunately I cant do this from my private computer so I have to wait until I am back at the institute... but many thanks already at this point.
Isa0984 is offline   Reply With Quote
Old 03-18-2015, 10:14 AM   #10
Isa0984
Member
 
Location: Germany

Join Date: Aug 2014
Posts: 15
Default

Hello, its long time ago, but still/again present for me... It was not possible for me to check the data again at the institute with samtools, but shoudn't I see the same (different read sizes) if I look with IGV to my data? That in fact gives me for all reads the same size of 51 bases, which means the campany didn't trimm the data before mapping... am I right? Thanks for your help! Isabelle
Isa0984 is offline   Reply With Quote
Old 03-18-2015, 11:07 AM   #11
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Yes, it sounds like they didn't trim them then. Scroll through IGV and see if there are any soft-clipped alignments (alignments that appear shorter but where the original sequence is 51). Using an aligner that does soft-clipping alleviates some of the issues surrounding adapter contamination and quality. If, however, they did end-to-end alignment (i.e., there are no soft-clipped alignments) on untrimmed data then I'd say they did a half-ass job.
dpryan is offline   Reply With Quote
Old 03-19-2015, 02:34 AM   #12
Isa0984
Member
 
Location: Germany

Join Date: Aug 2014
Posts: 15
Default

Hey, thanks for the fast replay. I found some shorter ones... they did it with the -q option of BWA.
When I asked them for the mapping parameters I got following answer:

n NUM max #diff (int) or missing prob under 0.02 err rate
t:4 (number of threads)
M:3 (mismatch penalty)
q: (quality threshold for read trimming down to 35bp 0)

I am not shure if I understand this 35bp thing, because I can find reads with a length less then 35bp (The 0 is maybe a typing error)?
Another question is, how can I get alignments like that (see figure)??? If you have n=0.02, shouldtn there at most 2 mismatches per 50 bp? Isabelle
Attached Images
File Type: png igv_snapshot.png (16.2 KB, 9 views)
Isa0984 is offline   Reply With Quote
Old 03-19-2015, 04:16 AM   #13
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Can't say I'm overly familiar with bwa aln, since most people use bwa mem these days.

The -n option has to have one of the more confusing descriptions I've seen. If it's an integer then the explanation is simple. I assume that it uses a poisson distribution with fractional -n, so a value of 0.02 with 50bp reads would correspond to a maximal edit distance of 3 (in R: qpois(0.98, 50*0.02)).

The -q option in bwa aln doesn't really specify a minimum read length. It specifies a value used when determining the trim location:

The -q value is INT and the quality at position i is q_i. So, this basically sums the penalties and finds the maximum value. The position with the maximum value is where trimming will occur (essentially, obviously if the penalty is <0 then no trimming should occur).
dpryan is offline   Reply With Quote
Old 03-19-2015, 05:40 AM   #14
Isa0984
Member
 
Location: Germany

Join Date: Aug 2014
Posts: 15
Default

Ok, I think I got the -q option, its just the information of the company, which is strange, maybe they mean a quality treshold of 35...
But the -n value is absolutely confusing... I was reading a lot of threads about this topic, but still. If in my case the maximal edit distance is 3, what does that mean??? Is there any relation to the allowed amount of mismatches?
Isa0984 is offline   Reply With Quote
Old 03-19-2015, 05:56 AM   #15
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

They are related, yes. "Edit distance" is a generalization of mismatches. If a read aligns with 3 mismatches then its edit distance is 3. However mismatches can't describe things like insertions or deletions. So if your read aligns with an insert of 2 bases then it has an edit distance of 2. If it has a single base mismatch and later a deletion of 3 bases then the edit distance is 4. The wikipedia article on edit distance is quite good. In short, "edit distance" is the minimum number of single character changes (insertion, deletion, or substitution) needed to convert one sequence to another.
dpryan is offline   Reply With Quote
Old 03-19-2015, 06:07 AM   #16
Isa0984
Member
 
Location: Germany

Join Date: Aug 2014
Posts: 15
Default

So it is the amount of operations (substitution, insertion, delition) that I have to do, to change my 50bp read to my reference. But why do I have then reads mapped to my reference with partly more than 10 mismatches?
Isa0984 is offline   Reply With Quote
Old 03-19-2015, 06:10 AM   #17
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

You'd have to go through the exact command that the company used and compare that to how things are processed inside bwa (by looking at the code rather than the documentation).
dpryan is offline   Reply With Quote
Old 03-19-2015, 06:18 AM   #18
Isa0984
Member
 
Location: Germany

Join Date: Aug 2014
Posts: 15
Default

That means other options in BWA can change the maximum of my 3 "mismatches" and allow more... wow that makes it even more complecated for me. They said they used default parameters with the only change in the options mentioned above.
Isa0984 is offline   Reply With Quote
Old 03-19-2015, 06:23 AM   #19
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,060
Default

You are not able to get the raw data and recreate the result yourself? You could be spinning your wheels for a long while trying to find an explanation for their results.

Note: Reading the beginning of the thread seems to indicate you have the raw data but are your results not matching the alignments the company did?

Last edited by GenoMax; 03-19-2015 at 06:26 AM.
GenoMax is offline   Reply With Quote
Old 03-19-2015, 06:29 AM   #20
Isa0984
Member
 
Location: Germany

Join Date: Aug 2014
Posts: 15
Default

I have the raw data, but I am not working at the institut anymore... and at home with my old little labtop, not really. Iam writing my thesis and there you look sometime at things deeper, then you should have done it before. But maybe I should simply take it like it is. Its just the urge to understand it in detail. Thanks a lot for your patience!
Isa0984 is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 04:37 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO