Mate pair sequencing - quality, duplication, throughput

[agupta29]

We recently got Illumina mate pair sequencing data for a human sample. I am trying to do some some Quality Check to see how good the data is.

I used bwa for alignment and the output bam file gave these numbers:
324061746 + 0 in total (QC-passed reads + QC-failed reads)
265547763 + 0 mapped (81.94%:-nan%)
239584614 + 0 properly paired (73.93%:-nan%)
244063409 + 0 with itself and mate mapped
21484354 + 0 singletons (6.63%:-nan%)
993167 + 0 with mate mapped to a different chr
812922 + 0 with mate mapped to a different chr (mapQ>=5)

I knew (from FastQC) that the data had high rates of duplication (>90%), so I marked and removed duplicates with picard to get these numbers from the dedupped bam file:

114026332 + 0 in total (QC-passed reads + QC-failed reads)
55512349 + 0 mapped (48.68%:-nan%)
114026332 + 0 paired in sequencing
57259162 + 0 read1
56767170 + 0 read2
49051344 + 0 properly paired (43.02%:-nan%)
50304349 + 0 with itself and mate mapped
5208000 + 0 singletons (4.57%:-nan%)
371423 + 0 with mate mapped to a different chr
246089 + 0 with mate mapped to a different chr (mapQ>=5)

As you can see, ~200 M reads are duplicates and consequently, removed from the alignment file.

My questions are:
1) Is it reasonable for human mate-pair libraries (insert size ~5kb) to have such high rates of duplication?
2) Does this reflect an average/good/bad mate pair sequencing run?
3) Any other suggestions regarding checking quality of mate pair sequencing data in general.

I have looked at other threads and this seemed like the only one somewhat relevant.

Thanks

[rhinoceros]

General related question: What is the largest insert size that can be used with Illumina sequencing if you want the pairs to actually overlap? Is it something like 200bp?

[luc]

Hi,
which MP library protocol did you use? One using linkers or not? Did you trim the linkers?
How exactly did you carry out the alignments?

[luc]

MiSeq does allow for 2x250 bp PE sequencing;
HiSeq does allow for 2x150 bp PE sequencing;

How much overlap do you want?

There is one program which claims to merge reads even without actual overlap - given enough coverage:

"COPE: An accurate k-mer based pair-end reads connection tool
to facilitate genome assembly"

General related question: What is the largest insert size that can be used with Illumina sequencing if you want the pairs to actually overlap? Is it something like 200bp?

[GenoMax]

General related question: What is the largest insert size that can be used with Illumina sequencing if you want the pairs to actually overlap? Is it something like 200bp?

Upwards of 400 bp. With MiSeq 2 x 250 bp runs you will get overlap in the middle. For HiSeq runs would be shorter since 2 x 150 bp is the max supported on the 2500 in rapid mode at the moment.

[agupta29]

The libraries were prepared by Illumina Mate Pair Library Preparation Kit v2 and sequenced using HiSeq2000 generating 2*100bp reads.

The alignments were done using bwa aln with default parameters. I did no pre-processing of data in this current analysis.

bwa aln -t 4  /BWAIndex/genome.fa L001_R1.fastq.gz > L001_1.sai
bwa aln -t 4  /BWAIndex/genome.fa L001_R2.fastq.gz > L001_2.sai
bwa sampe /BWAIndex/genome.fa 1.sai 2.sai 1.fastq.gz 2.fastq.gz > L001.sam
samtools view -bS L001.sam > L001.bam
samtools sort L001.bam L001_sorted

The duplicates were marked and removed, as I previously mentioned, using picard MarkDuplicates tool.

Hi,
which MP library protocol did you use? One using linkers or not? Did you trim the linkers?
How exactly did you carry out the alignments?

[luc]

Hi,

I believe the "old" bwa has problems with mapping mate pairs:

http://seqanswers.com/forums/showthread.php?t=5085

BWA-mem, Bowtie, Novalign are better suited to my knowledge.

Similarly the old Illumina MP kit produces lots of artifacts and PCR duplicates; the Nextera based mate-pair kit is more efficient.

What was the fragment size of your library? People used to trim the MP reads to the first 38 bases each, in part to avoid mapping chimeric sequences.

[agupta29]

BWA and BWA mem
I an using reverse-complemented reads. That brings the --rf orientation mate pair reads into --fr paired end orientation, which can then be used for bwa alignment. I am running bwa-mem as I write this so that I have a direct comparison between bwa/ bwa mem

Library prep:
The samples were sent to a commercial vendor. So, I did nothing with the library or sequencing steps. We just got the data from the vendor. But, I will keep that in mind for the future.

[luc]

Hi agupta,

you could look at the distances of PE alignments to get an idea if their sequencing libraries had short inserts and thus were more likely to include chimeric reads.

[diptarka]

I have carried out de novo assembly of an organism of interest using velvet optimiser. I have contig files. How can i now, predict genes from the sequences? I have also carried out comparsion of the assembly with a reference using abacas and tried to view it usin artemis comparsion viewer. But, i am unable to understand the contig ordering part. How does one do that? Secondly, what is the way to predict genes from the contigs?

[GenoMax]

I have carried out de novo assembly of an organism of interest using velvet optimiser. I have contig files. How can i now, predict genes from the sequences? I have also carried out comparsion of the assembly with a reference using abacas and tried to view it usin artemis comparsion viewer. But, i am unable to understand the contig ordering part. How does one do that? Secondly, what is the way to predict genes from the contigs?

You should create a separate thread for this question since:
a) your questions are not related to current thread
b) people who could potentially answer your question will not notice it if it is embedded here.

New threads can be started by doing following:

Seqanswers.com --> Forums (Navigation Menu on Left) --> Select an appropriate forum (e.g. Bioinformatics) --> "New Thread" button at the top left corner on next page.

That said: What type of organism is this (prokaryote/eukaryote)?