I’ve talked
about FastQC in a
previous post but now I’ve actually use it to analyze our data. We want to make sure the data
we have is of sufficient quality before we do anything with it, so it’s
important to do this analysis early.
Note: I put all relevant documents on the Google Drive
in the folder Scott>Quality Analysis
Here’s what
I did:
1)
Run FastQC on all .fastq files (96 x 2 reads per
sample = 192 files)
2)
Run Sambamba
to align files to reference genome (to be discussed in future post)
3)
Run Sambamba
to extract alignment statistics
4)
Made a script that runs in the terminal to
extract important statistics from the output of 1) and 3) (uploaded to the
drive as fastqc_extraction.sh)
The output
of my script is called fastq_results.xlsx on the drive. On the first worksheet
of this file I have listed useful values and whether or not a sample passes or
fails the standard FastQC tests. The second worksheet is just a list of every
sequence FastQC has labelled as “overrepresented”.
The first
thing I did was compare some numbers between the old and new RNA seq datasets
in R. The script that makes this plot is called quality_analysis_script.R (and requires
the script multiplot.R)
Number of reads (i.e. number of DNA
fragments sequenced in each sample)
The graph
suggests that the older samples tend to have more reads. Older samples have, on
average, 1.5 million more reads. A t-test says that there is a significant
difference with a p-value of 1e-08. Although 1.5 million sounds like a lot,
keep in mind that some of the older samples differ by close to 6 million reads.
The top 4
files for reads are:
toxx_M0_A.R1.fq.gz
– 10651839
toxx_M0_A.R2.fq.gz
– 10651839
antx_M0_A.R1.fq.gz
– 11311458
antx_M0_A.R2.fq.gz
– 11311458
The bottom 4
files are:
antx_M2_C.R1.fq.gz – 1781164
antx_M2_C.R2.fq.gz – 1781164
antx_M2_E.R1.fq.gz
– 4035388
antx_M2_E.R2.fq.gz
– 4035388
Number of duplicates (i.e. reads that do not
have a unique sequence in the sample)
The graph
suggests that the newer samples tend to have more duplicated reads. Newer
samples have, on average, 4% more duplicates. A t-test says that there is a
significant difference with a p-value of 4e-11.
Duplicate
reads are caused by sequencing two PCR products derived from the same fragment
(technical duplicates) or by sequencing the same region of two distinct RNA
molecules (biological duplicates). Technical duplicates are unwanted whereas
biological duplicates are pretty much the basis of differential expression
analysis. It is impossible to discriminate the two types of duplicates with the data we have.
So it looks
like most samples fall above 80% duplicates – which sounds like a lot (only 20%
of reads are actually unique). However, remember that only one of the two paired
end reads is actually being analyzed here so all reads that begin at the same nucleotide
are considered duplicates. Because we’re doing RNA seq here, most reads come
from the transcriptome (a subset of the whole genome) and in some samples,
there are 3x more reads than nucleotides in the genome; extremely high duplication
(80-90%) is reasonable.
The top 4
files for duplicated % are:
kw20_M3_C.R1.fq.gz
– 89.1
toxx_M1_D.R1.fq.gz
– 88.6
kw20_M3_C.R2.fq.gz
– 88.0
toxx_M1_D.R2.fq.gz
– 87.9
The bottom 4
files are:
kw20_M0_B.R2.fq.gz
– 75.8
sxyx_M1_B.R2.fq.gz
– 75.5
antx_M2_C.R1.fq.gz – 60.4
antx_M2_C.R2.fq.gz – 58.7
Note the
bottom two. They have a lot of contaminated sequences, so it’s not surprising that
almost half of sequences are unique.
GC Percent (i.e. % of nucleotides that are G
or C)
There is on
average 1% more GC in the old data than the new (p=2e-08). H. influenzae has a GC % of 38% and Rosie has claimed that genes
tend to be higher. It seems that the older dataset is much more variable
(except for 4, 2 of which have been shown to be contaminated so far)
The top 4
files for GC % are:
kw20_M3_C.R1.fq.gz
– 46
kw20_M3_C.R2.fq.gz
– 45
antx_M2_C.R1.fq. gz – 45
antx_M2_C.R2.fq.gz – 45
The bottom 4
files are:
Most files
are lowest at 40%
% Alignment (i.e. % of reads that can be
mapped back to the H. influenzae genome)
Note that
for antx_M2_C, only 20% reads aligned to the reference genome. This sample was not included in the graph. On
average, reads from the new samples get aligned better than reads from the old
sample by 0.5%. This is significant (p=1e-10) once the poorly aligned sample is
removed.
This means
that the newer samples could be purer or the quality of sequencing was higher
in the new sample (base calling errors make it harder to align back to a
reference sequence) or maybe something else I haven’t thought of.
The top 4
files for GC % are:
sxy1_B2_G.R1.fq.gz
– 99.42
sxy1_B2_G.R2.fq.gz
– 99.42
sxyx_M0_D.R1.fq.gz
– 99.43
sxyx_M0_D.R2.fq.gz
– 99.43
The bottom 4
files are:
antx_M2_C.R1.fq.gz – 20.1
antx_M2_C.R2.fq.gz – 20.1
toxx_M2_A.R1.fq.gz
– 97.1
toxx_M2_A.R2.fq.gz
– 97.1
Overall:
-The old
data and the new data are different.
Statistics can tell me that, however I don't see any major problem. Sample antx_M2_C is what I would expect a problem would look
like. When samples are prepared by different people on different days, there is
bound to be technical variation.
- Sample antx_M2_C is bad. I feel it may cause problems to keep
in in any future analysis steps (but we can try anyway)
- Looking at
the graphs, my prediction is that the newer data is either less contaminated or contains
more PCR duplicates, because the graphs suggest that there is more H. influenzae DNA in the newer sample
(better alignment, GC content has low
variability)
So far, this
analysis suggests that the new data is probably as good as (if not better) than the old
data. In the next post I will look at the FastQC test results.
