constitution of read and gene region (IGV)

Question

I work with fastq files containing NGS reads for some human DNA regions. The reference genome is hg19. I had two fastq files (pair-ended). I generated alignment BAM files. I used "bwa" and samtools to find a possible target gene region (chr7:55,242,376-55,242,574,). This corresponds to a region of the EGFR gene.

Here is a screenshot of the gene region

And here is a screenshot of the primer region

I also have a list of primers and reverse primers:

FORWARD
1. TTGCCAGTTAACGTCTTCCTTCTCTCTCTG
2. CCCTTGTCTCTGTGTTCTTGTCCCCCCCA
3. TGATCTGTCCCTCACAGCAGGGTCTTCTCT
4. CACACTGACGTGCCTCTCCCTCCCTCCA

REVERSE
1. GAGAAAAGGTGGGCCTGAGGTTCAGAGCCA
2. CCCCACCAGACCATGAGAGGCCCTGCGGCC
3. TGACCTAAAGCCACCTCCTTA
4. CCGTATCTCCCTTCCCTGATTA

And I have some adaptors

ADAPTORS 1
AAGACTCGGCAGCATCTCCA

ADAPTORS 2
GCGATCGTCACTGTTCTCCA

I have to answer to the following questions:

1) What is the constitution of each read (adaptor+primer+amplified region)? The forward primer is obvious and corresponds to the first primer of the list

TTGCCAGTTAACGTCTTCCTTCTCTCTCTG

So why do we have the three other forward primers? I don't understand. And what about the correct reverse primer?

It seems to be the complementary sequence of the end region

GAGAAAAGGTGGGCCTGAGGTTCAGAGCCA

Finally which adaptor is used?

Is [ADAPTOR1 - PRIMER1 - AMPLIFIED REGION] the correct answer? Are there other possibilities?

2) % of reads that map the human genome? what about the unmapped?

Can you please give me some clues about how to answer this question.

I just need some clues I want to do it by myself.

Thank you very much for your help.

Answer 1

You have asked quite a few questions here. I will try to answer a subset of them.

To answer your question, I started by going to the UCSC genome browser, and selected the BLAT from the tools dropdown menu. Then passing your primer sequences as queries, we can clearly see where they map to the reference genome.

Looking at the figure above, I can guess that your data comes from a targeted sequencing experiment (exome, or gene panel, etc); where you have primers around each exon of EGFR on two strands, on 5' end of exon sequence on each strand (think about DNA polymerase acting in PCR to amplify the number of DNA molecules before sequencing).

Regarding the adapter, I am not sure which read you are referring to. Can you clarify please?

About percentage of mapped vs unmapped reads, you can use bamtools.

Example Usage:

/user/me/src/bamtools/bin/bamtools-2.3.0 stats -insert -in my-sequence-file.bam

Example output:

Total reads: 103277668

Mapped reads: 90088436 (87.2293%)

Forward strand: 58136735 (56.2917%)

Reverse strand: 45140933 (43.7083%)

Failed QC: 6529806 (6.32257%)

Duplicates: 0 (0%)

Paired-end reads: 103277668 (100%)

'Proper-pairs': 87439672 (84.6646%)

Both pairs mapped: 87910438 (85.1205%)

Read 1: 51638834

Read 2: 51638834

Singletons: 2177998 (2.10888%)

Average insert size (absolute value): 6317.39

Median insert size (absolute value): 301