Exome sequencing alignment

[Jackken]

Hi,

I used bowtie to align exome sequencing (Illumina GA), and this is what I got:


# reads processed: 37205349
# reads with at least one reported alignment: 141065 (0.38%)
# reads that failed to align: 37064284 (99.62%)
Reported 141065 alignments to 1 output stream(s)

I am wondering what went wrong.

By the way, here is the .fastq file looks like,

@SRR350953.3 MENDEL_0047_FC62MN8AAXX:1:1:1488:946 length=152TTTTTTTT
NTCCCATTATCTCAAGCAGCCATATGTTTCTCATTCACTTGATACACTGTTTCTTTTCAACCCCCACATCCTCACCGTGCTCAA
ACAAAGAAACAGGTGGTGAGGATGTGGGGGTTGAAAAGAAACAGTGTATCAAGTGAATGAGAAACATA########B@<7:EFE
+SRR350953.3 MENDEL_0047_FC62MN8AAXX:1:1:1488:946 length=152EEB@@>>D
############################################################################FAD=FCCC
CCFDD;E??@@FD?BD>F?FB=BFAEFDGGDGGDG@BB=5=/;?8B=DDDGDGBDGDACAE@B@G88;CTCCTTCCAGGACCCA


Thanks!

[Jackken]

Exome sequencing alignment

Can anybody help?

Thanks!

[TonyBrooks]

"#" is a Q-score of 2 (>60% error). There's your problem.

Actually, after that the error rate improves with still plenty of sequence left in the read. You can try trimming your reads using fastq-trimmer and re-aligning.

[Jackken]

Thanks! I will give it a try.

[Jackken]

fastx_trimmer gave me an error (see below), is there anyway to make it work?

Thanks!

fastx_trimmer: Invalid quality score value (char '#' ord 35 quality value -29) on line 4.

@SRR350953.1 MENDEL_0047_FC62MN8AAXX:1:1:1206:930 length=152
NTGATTTAGCTGCATAGTTTTCTTCTTTTTAATCCATAATGTATACATTTTAGACTTTGTATTTTAACTGCTGACATTCC
AGTCTAAGTCGGAAGCCACATCTTCTAAACCAAATGTCTCTTCATCCCTTATGTCAGGAACCTATTTTTTTT
+SRR350953.1 MENDEL_0047_FC62MN8AAXX:1:1:1206:930 length=152
############################################################################B@<7
:EFEEBF?8B?3=;@9GGGG?;:C7CBABA=DG><GGB>DGE>3<EADGEC=DDB8GGD3<CE-EEB@@>>D

[quique_vzquez]

Are you sure your data isn't paired end? When I've got that large reads from Illumina always are paired end.
If it is pair, you must separate your file into two files before align.

[Jackken]

Many thanks!

I am now trying to use "grep" to separate the original file into two.

grep -A 1 "\.1 " originalfile.fastq > newfile_1.fastq
grep -A 1 "\.2 " originalfile.fastq > newfile_2.fastq

[quique_vzquez]

I separate them from the original .sra file with:
fastq-dump --split-3 originalFile.sra

[Jackken]

Thanks a lot!

I will try it.

[kmcarr]

fastx_trimmer gave me an error (see below), is there anyway to make it work?

Thanks!

fastx_trimmer: Invalid quality score value (char '#' ord 35 quality value -29) on line 4.

@SRR350953.1 MENDEL_0047_FC62MN8AAXX:1:1:1206:930 length=152
NTGATTTAGCTGCATAGTTTTCTTCTTTTTAATCCATAATGTATACATTTTAGACTTTGTATTTTAACTGCTGACATTCC
AGTCTAAGTCGGAAGCCACATCTTCTAAACCAAATGTCTCTTCATCCCTTATGTCAGGAACCTATTTTTTTT
+SRR350953.1 MENDEL_0047_FC62MN8AAXX:1:1:1206:930 length=152
############################################################################B@<7
:EFEEBF?8B?3=;@9GGGG?;:C7CBABA=DG><GGB>DGE>3<EADGEC=DDB8GGD3<CE-EEB@@>>D

The FASTX toolkit still assumes by default that all FASTQ files use the original Solexa Phred+64 encoding for their quality scores. Your file uses the (now standard) Phred+33 encoding. You have to explicitly tell fastx_trimmer that your file is Phred+33 by adding the parameter "-Q33" to your command line.

[Jackken]

Thanks, kmcarr.

I think I didn't realize that it's paired end. So quique_vzquez is right. And I am separating the original .fastq file into two. I think it's working now.

quique_vzquez, thanks a lot!