Umbalanced content of complementary bases in DNA-seq paired-ends

[ratope]

Hi Community:

I am doing DNA-seq (paired-ends) of a eudicot plant. FastaQC shows "per base sequence content" graphics (one per each "paired-end" fastaq file) that puzzle me. As you see, %A and %T are very different and %C and %G are also different. On the other side, the two graphics are in a sense "complementary" i.e. A1+A2 ~= T1+T2 and C1+C2 ~= G1+G2.

Additionally, I have counted the nucleotides in both entire "paired-end" files, using only the reads that aligned after BWA. Please, compare it with the same figures in the complete genome:

1_1 reads (fastaq file):
A 1454788765 32.45%
T 1007930585 22.48%
C 875193373 19.52%
G 1109913625 24.76%

1_2 reads (fastaq file):
A 1051644195 23.46%
T 1422158391 31.72%
C 1075609393 23.99%
G 896848970 20.01%

Complete genome (deposited fasta file):
A 153891618 31.65%
T 153909273 31.65%
C 81246927 16.71%
G 81219831 16.70%

Other previous DNA-seq paired-end experiment showed A=T and C=G in each single file (see attached image num. 3). It seems to me the logical result, as I presume each subset of files contain a not-biased set of reads.

Have you seen this kind of results before? Have you an explanation of both: a) the umbalance between the complementary bases and b) the complementarity of the results of both paired-end files?

Thank you in advance!
Rafael

[kmcarr]

How was the sequencing library prepared? Your numbers suggest that the library was prepared in some sort of strand specific manner.

[ratope]

Thank you for your interest, kmcarr.

As far as I know, no strand specific protocol was used by the institution that made the sequencing. Additional data: 2 X 150 pair-end sequencing on a GA-IIx. 1 sample in 2 lanes (channels). Insert size: 400 bp.

[fkrueger]

To me it looks like your sequencing captured (for whatever reason) reads with a higher proportion of C and Gs, which are poor in As. While I don't know the reason for this shift in base composition (possibly over-amplified CT-rich repeat regions?), it seems to make perfect sense that the you are seeing an inverse pattern on the paired-end side. Let's say you sequence a fragment which has a C content of ~25% and a T content of 35%, then the reverse complemented sequence (what read 2 sequences) would a similarly high G or A content, respectively, and this seems to be quite well the case.

The base composition in the raw sequence data may very well be completely different from the composition you counted for the mapped data. This can be due to contamination of the sequencing sample by DNA of another species, adapters or sequences not yet in the assembly to name but a few. Did you see any hints of contamination in the %GC or k-mer plots?

[kmcarr]

To me it looks like your sequencing captured (for whatever reason) reads with a higher proportion of C and Gs, which are poor in As. While I don't know the reason for this shift in base composition (possibly over-amplified CT-rich repeat regions?), it seems to make perfect sense that the you are seeing an inverse pattern on the paired-end side. Let's say you sequence a fragment which has a C content of ~25% and a T content of 35%, then the reverse complemented sequence (what read 2 sequences) would a similarly high G or A content, respectively, and this seems to be quite well the case.

But the standard protocols for preparing Illumina libraries are not strand specific. That is to say it is a random selection which strand of a sheared and polished fragment will be read 1. Take this hypothetical as an extreme example, a collection of fragments which are simply the repeat AACAAC ad infinitum. The complementary strand would be TTGTTG. Since it is equally probable that either strand would end up as read 1 then you still should see balanced (A==T, G==C) base composition.

[ratope]

Thank you for your suggestions, fkrueger and kmcarr!

Yes, the GC content is considerably increased (~9%) respect to the theorical for this species. There are some overrepresented sequences (long N..N sequence and a primer of the paired-end) but in a low proportion.

I have went forward and in fact there is a percentage of 51% of reads mapping to the reference. Doing blastn on a small sample of non-mapping reads I find that most of them (~92%) have hits in both contigs (the majority) and chloroplastic/mitochondrial DNA of my species. My impression with these data is that there is no substantial contamination, but perhaps I am wrong...