SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
BFAST for SOLiD paired end reads epigen Bioinformatics 31 09-03-2011 05:20 AM
BFAST mapping paired end reads. tanghz Bioinformatics 10 04-29-2011 06:29 AM
How to map SOLiD paired end reads by Bfast beliefbio Bioinformatics 1 12-29-2010 12:55 AM
BFAST input format for paired end reads lindseyjane Bioinformatics 5 12-16-2009 07:21 AM

Reply
 
Thread Tools
Old 01-04-2012, 02:58 AM   #1
gavin.oliver
Senior Member
 
Location: uk

Join Date: Jan 2010
Posts: 110
Default Using Bfast to align paired end Illumina reads

Hi,

I am attempting to align paired end Illumina reads to the human genome with Bfast.

The problem is that most aligners expect paired end reads in 2 fastq files while Bfast expects 1 file. Having searched this and other forums, the most commonly suggested/accepted solution appears to be conversion of the 2 files to 1 file using the ill2fastq script supplied with Bfast.

I have now converted the files and run the alignment however there appear to be problems with the results. I have aligned the same sequences with Stampy and BWA (separately) and I get completely different results. The output of samtools flagstat is shown below.

Does anyone have any suggestions on how best to utilise Bfast under such circumstances?

BWA
200000 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
191942 + 0 mapped (95.97%:nan%)
200000 + 0 paired in sequencing
100000 + 0 read1
100000 + 0 read2
190714 + 0 properly paired (95.36%:nan%)
191364 + 0 with itself and mate mapped
578 + 0 singletons (0.29%:nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

Stampy
200000 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
199376 + 0 mapped (99.69%:nan%)
200000 + 0 paired in sequencing
100000 + 0 read1
100000 + 0 read2
197774 + 0 properly paired (98.89%:nan%)
198768 + 0 with itself and mate mapped
608 + 0 singletons (0.30%:nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

BFAST
125396 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
123238 + 0 mapped (98.28%:nan%)
50792 + 0 paired in sequencing
25396 + 0 read1
25396 + 0 read2
0 + 0 properly paired (0.00%:nan%)
49516 + 0 with itself and mate mapped
326 + 0 singletons (0.64%:nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
gavin.oliver is offline   Reply With Quote
Old 01-04-2012, 03:55 PM   #2
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

If you have two files, then you need to create one file that is interleaved so that the second read of a pair is directly after the first. Then specify the pairing orientation on the command line.
nilshomer is offline   Reply With Quote
Old 01-05-2012, 04:15 AM   #3
gavin.oliver
Senior Member
 
Location: uk

Join Date: Jan 2010
Posts: 110
Default

Thanks Nils - I'll write something to interleave the two files.
gavin.oliver is offline   Reply With Quote
Old 01-06-2012, 12:30 AM   #4
gavin.oliver
Senior Member
 
Location: uk

Join Date: Jan 2010
Posts: 110
Default

Nils,

I interleaved the two files and have run them using Bfast. However, whenever I reach the postprocess stage, I am not getting successful pairing of the reads regardless of what I specify on the command line with the -P/-S or -Y options.

Here is a small selection of the reads:

Code:
@chr17:95861:F:237/1
GTTGGCGAACACATCCATGTGCCGGGAGGATGGTGCACCCCAACTCCACAAGGACCCTTCCAGACCGCGGCCGCTCCAGCTCTCAAAGCC
+
@<@?DDDD>FF:DAF9FFFCAGF<F3AAFD>2ACEF?CFC@?;FB:?@?;D@>86';EE;AE376?########################
GCCCCAGACTCCACAGGTTAAGGGCTCGCATCTCTTGAACAGGGATCTTGATTGCCCCGCGACCTACTGACAATCTGAATTCTGGGCGCT
+
+:+=+ADDFCFBF3C:BE+@@;1):C?###############################################################
@chr17:42735:F:247/1
TTAAAACTGGATATCACCCAGTGTTGGCAAGGTACAGGAAAATGGGAACTATCATATACCACAGGGGCTGGAAGAGCATAAACTGGTTTA
+
CCCFFFFFHHHHHJJJJJJJJIJJJJJJJJJIJGIJJJJJJJIJJJJIJJJJJJIJJJJJIJIJIJJJIJIIJJJJIIJJIJHGFHHFFF
@chr17:42735:F:247/2
TTAAGTGACTTCATTTTTAATTACTATATGGGATTCTATCTTTCCAGTGTATCATGATTTATTTGACCTATTGCTGAATGTTGGAGGTTT
+
@CCFFFDDHHHHGJJJJJHIJJJJIGJJJJCHJJJIHIJIJJJJJIIIIJIIIGIIJJJHHHHIIGEE@HC=?BE>;>CEEC>A=ACC;A
@chr17:85755:R:-161/1
AAAATAATAAACCAGTCATTAGAACCATAAACCTGTACTGTTTTTGACAATGTAATGCACTGCCCTGTAAAGCACTACAATAAAGACGTT
+
1:1+4??=DCFC>+A)++3+3CACG4??:*::?DCDEAGB1?DB?<?B<B@DEHFE=)8B8=77=@((6).=CC7=;@>D@3)7;>B###
@chr17:85755:R:-161/2
AATTACATGTCATATGAATGAATACTTGACGTCAGCAGGACTGCGTTTTGGTGGTGAACTTGGTTCTAGGTAGAAACAAAGAATGGAGAA
+
@@@DAD:DFDFFDF9A4CB<:A?:<FC)ACF6?<CFFF;GABAB93BDGCE<DB?*8=@=B@=FFGIF<1?@>?BB@;AC>;(5-;;>A>
@chr17:76132:F:140/1
AATGTAAAACACCATAAAAATTAATCTTAAGGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGGGAGGCCGAGGTGGGCGGA
+
CCCFFFFFHHHHHJJJJIJJIJJJJJJJJJJJJJJJJJJJJIJIJJJJIJJJJGIJJIJHIGJJJHGIJJJJJJJIHGIIFIHHHHHHFF
@chr17:76132:F:140/2
ACGCCATTCTCCTGCCTCGGCCCCCCAAGTAGCTGGGACTACAGGCGTCCACCACCACGCCCGGCTAATTTTTTGTATTTTTAGTAGAGA
+
@C@FFFFFFHHHGJJJJJJJJJJJJJJIJJJJEHIIJJJJJJJJJIJDGHGGIIGHIIIJJJJJJJJGIIJJJIHGHHHHFFFFFDEEEE
@chr17:72239:F:208/1
CAGACTGGCATGCAATGGTGCGATCTCGGCTCACTGCAACCTCTGCCTCCCAAGTTCAAGTGATTCTCCTCCCTCAGCCTCCCAGGTAGT
+
CCCFFFFFHGHHHJJJJJJIJJGIGEHGIIJJHIJJJJJJJJJJJGJJJJIJIJJIIJIJEGCGHIIIHHDHFHHHFFFFFFFDCEEECC
@chr17:72239:F:208/2
AGAAAAAGTAAGAGACACCTATAGATCAGAAGACACTTGGGGCTGGGCATGGTGGCTCACACCTGTAATCCCAGCACTTTGGGAGGCCAA
+
@CCFFDDFHHHHHJJJFHIJJGJIJJJJJEGGIIJJJJJGIIJJJJGEIIFHGHIFHIIJJJJJEIJIIJIIGHHFFEEHFFFFFFEFCC
@chr17:124386:R:-215/1
AAGGGAAAAGACCCAAAGGGTTGGAAGCAATATGTGAAAAAATACAGAATTTATATTGTCTAATTACAAAAAGCAACTTCTAGAACCTTT
+
CC@FFFFFHHHHGJJJJJJJJJJJJJJJJJJJJJJIJJJJJJJIJHJJIJHIJJJJJJJJJJJIIIJJIIGIJIJIJGJJJIHHHHHHFF
@chr17:124386:R:-215/2
AGCCAGGCATGGTGGTGCATGCGTGTAATCCCAGCTACTCGGGAGCTGAGGCAGGAGAATCGCTTGAACCCACGAGTCAGAGGTTGCGGT
+
CCCFFFFFGHHHHJJJJJJHIJJJJJJIJIIIJJJIJJJJGIJIJJJJJJIJJIJJJJJJIIJJJIJJJHHHHHHHHFFFFFFFDFEEED
@chr17:43455:F:112/1
TTTTTGAGACAGGGTCTCACCATATCACCCAGGATGGAGTGCAGTGGCACCATCATGGCTCACCACGGCCTCAACTTGCTGGGATCAAGC
+
CCCFFFFFHHHHHJJIJGIJJJIJGIJJJJJJJJJJIIIJJJJIIJJJJJIJFGHGIJJJJEIIGHHHGHFFFFFEEEEEEEDDDDDDDC
@chr17:43455:F:112/2
ACGCTACTGCACTCTGTTCTAGGCAACCCCTGTCTGGGAAAAAAAAAAAAATTAGTGAGGCTTAGTGGTGCACACCTGTAGTCTCAGCTA
+
CCCFFFFFHHHHHJJJJIEIJJJJHGIJJIIEHGIIJJIIJ@FHIJJJIJGHIJIHHBGHGI=FIEDEGGIJHEHG=>AE?EHCEFDFE;
@chr17:89085:R:-274/1
TTTAAAGCACTACCCAGAGTTATTCAAAGCCAGGCAGGAGAACTGCAGGCATCAGAATGCCCTGGGGACGGGTCCAAAATGCAGAATCCT
+
@@BFFFFDHDHHHFGIJIGJJGJIJIHIEGJJJ)?@FHHJJI;FDHIIHH=>DFEDECEEECDDECCD55<@B>CDCBDD@CDCCDACCC
@chr17:89085:R:-274/2
GCTTGGTAAGGAATATTAGGTGAACACAGGTGCCTATGTGGGTTCCACCGTTTTCCATAAATGTTAATTTTCGGCACTGAGGACGGACAG
+
CCCFFFFFHHHHGIJJJJJIIIIJIJHF>CCE@3;(5>;?B?3:,5,3?((+3@ADCD::8:::4>@CD#####################
Is there anything obviously problematic to be seen?
gavin.oliver is offline   Reply With Quote
Old 01-06-2012, 04:30 PM   #5
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

That doesn't look like a valid FASTQ file (where's the read name for the second read). Also, the read names should be the same for multi-end reads. Please read the manual on the topic of the proper input format in the manual accompanying BFAST. There should be an example paired end read FASTQ (etc.).
nilshomer is offline   Reply With Quote
Old 01-08-2012, 11:52 PM   #6
gavin.oliver
Senior Member
 
Location: uk

Join Date: Jan 2010
Posts: 110
Default

Nils, RE the second read name - it must have been a copy/paste error as the file actually does contain it:

Code:
@chr17:95861:F:237/1
GTTGGCGAACACATCCATGTGCCGGGAGGATGGTGCACCCCAACTCCACAAGGACCCTTCCAGACCGCGGCCGCTCCAGCTCTCAAAGCC
+
@<@?DDDD>FF:DAF9FFFCAGF<F3AAFD>2ACEF?CFC@?;FB:?@?;D@>86';EE;AE376?########################
@chr17:95861:F:237/2
GCCCCAGACTCCACAGGTTAAGGGCTCGCATCTCTTGAACAGGGATCTTGATTGCCCCGCGACCTACTGACAATCTGAATTCTGGGCGCT
+
+:+=+ADDFCFBF3C:BE+@@;1):C?###############################################################
I'll have a look at the FASTQ format in the Bfast manual and reattempt. I hadn't scrutinised the format too deeply as the same sequences have already been successfully aligned with a few other programs.

EDIT: I have now rerun and the alignments look good. The problem seems to have been the presence of the /1 and /2 suffixes in the read names.

Last edited by gavin.oliver; 01-09-2012 at 12:56 AM.
gavin.oliver is offline   Reply With Quote
Old 01-09-2012, 03:23 PM   #7
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

Great, I am glad it is working for you!
nilshomer is offline   Reply With Quote
Old 01-11-2012, 01:28 AM   #8
david.tamborero
Member
 
Location: spain

Join Date: Feb 2011
Posts: 60
Default

Hi,

I have actually the same problem. I have each pair end of the reads in separate files, but when they look as the following:

--- pair_1 file:

Code:
@HWUSI-EAS1692_0001:1:1:1050:4451#0/1
CAGATTCACANTCCTGAATATCATGTTTTCTTTCCAAGGNATGACATAACGTCTTGGGATCATCCCTTGCTTTAATGAAAATCGTGGCAAATGAA
+HWUSI-EAS1692_0001:1:1:1050:4451#0/1
Ybaac][T^YB[ZZ[SKVZT`bcYbccaccaaa_cZZ[ZB[Z[T_c`cYcc\bcccc^T\a`TcccbL\ac\^a\Ybb`^bY]bb_BBBBB
-- pair_2 file:

Code:
@HWUSI-EAS1692_0001:1:1:1050:4451#0/2
CATGATAATGCACTCCATCTCATTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGCACTAAAAAGCGGACCTTGGTGTGAAAACATAACACACAC
+HWUSI-EAS1692_0001:1:1:1050:4451#0/2
M_M^ZM\YL]U^L\^VQJIU\a__\``c\cW_aaaaa_R[_\_`W][__BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
I use the ill2fastq.pl script contained in the bfast package to merge them and everything works fine. However, when the files are like:

-- pair_1 file:

Code:
@ILLUMINA-A16956_0001:4:1:1099:4197#0/1
CTGTTTTCAAAATAATATNCTTCTTGATTCCTTAGCTCTGTGCCTCAGAC
+
CCACCCACCCCCCC@<BA#@;:<<<CA=CCA=BBB@@C@ACCCCBCCCBC
-- pair_2 file:

Quote:
@ILLUMINA-A16956_0001:4:1:1099:4197#0/2
NNTATTTTCTAAAGTGAGGCGGCANGNANATNNTNGTAGGTANCNNNGNN
+
###########################################################
which seems more 'fastq - like' (with Phred+33 quality and '+' in the third line), then I have to merge them by my own script? I mean, I can't not merge them by using some argument of the bfast commands used for alignement?

Thank you very much for your feedback,

best
David
david.tamborero is offline   Reply With Quote
Old 01-11-2012, 05:11 PM   #9
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

The read names should be identical; they are not.
nilshomer is offline   Reply With Quote
Old 01-12-2012, 08:49 AM   #10
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

The de novo assembler velvet also comes with a little perl script to interleveave two separate fastqs into one.
swbarnes2 is offline   Reply With Quote
Old 01-14-2012, 06:12 AM   #11
david.tamborero
Member
 
Location: spain

Join Date: Feb 2011
Posts: 60
Default

Nils, sorry for that, but I would say that the reads name are the same. In the first example, they are 'illumina-like'. In the second example, they are 'fastq-like'.

So my guess is that in the first case they must be interleaved + converted, I did it by using the ill2fastq script contained in the bfast package, and everything went fine.

In the second example, they have to be interleaved (but not converted, since the quality is already +33).

So, sorry if I am missing something obvious, but my questions are the following:

- I noticed that the ill2fastq script also complements + reverts the sequence of the second pair. Is it necessary?

- is there any way to input the two fastq files to the bfast command(s), without having to interleave them before?

thanks,
david
david.tamborero is offline   Reply With Quote
Old 01-14-2012, 06:22 AM   #12
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

No, BFAST requires the reads to be interleaved with the read names exactly the same. Supporting platform specific vendors input formats becomes difficult (look at all of these conversion scripts!), and so a generic format is better.

I haven't looked at ILMN data in a while, so I am not sure if the ill2fastq scripts needs updating. I would be immensely grateful for a patch.
nilshomer is offline   Reply With Quote
Old 01-14-2012, 06:43 AM   #13
david.tamborero
Member
 
Location: spain

Join Date: Feb 2011
Posts: 60
Default

Yes, you're right.
Last question: when interleaving, the second end must be reversed and complemented?

I am trying to figure out whether bfast requires all the pairs to be listed in order of 5' to 3' and on the same strand. I am using bfast+bwa-0.7.0a.
david.tamborero is offline   Reply With Quote
Old 01-14-2012, 06:47 AM   #14
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

Prior versions required that, but now you can specify their orientation during the "postprocess" step.
nilshomer is offline   Reply With Quote
Old 01-14-2012, 06:51 AM   #15
david.tamborero
Member
 
Location: spain

Join Date: Feb 2011
Posts: 60
Default

ok, I will download the last version to simplify the interleave step.

Many thanks for your help, Nils.
david.tamborero is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 03:59 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO