SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
samtools mpileup calls way too less SNPs TuA Bioinformatics 17 03-01-2018 04:17 PM
trouble repeating a BWA alignment command efoss Bioinformatics 0 08-09-2011 01:51 PM
Bowtie vs BWA: only 50% overlapping SNPs a11msp Bioinformatics 4 10-14-2010 03:22 AM
Samtools pileup SNPs questions Melissa Bioinformatics 2 03-24-2010 10:15 AM
bwa - samtools - SNPs Mr Mutundes SOLiD 0 12-10-2009 03:11 AM

Reply
 
Thread Tools
Old 12-12-2010, 10:25 AM   #1
mindlessbrain
Junior Member
 
Location: qatar

Join Date: Dec 2010
Posts: 7
Default Trouble getting SNPS from bwa/samtools

Hello all,

I have a segment from a yeast, and a yeast genome to compare it to. I want to get a list of SNPS, and eventually find a known SNP in the yeast sequence.

When I run this through bwa and samtools it seems to work fine up until its time for me to get the SNP list. When I use the pileup command I get empty files as output.

Here's what I'm running.

Code:
#reference is yeast.nt
#query is yeast.fasta

bwa index /home/jill/bwa-0.5.8c/Project/yeast.fasta 

bwa index /home/jill/bwa-0.5.8c/Project/yeast.nt

bwa aln /home/jill/bwa-0.5.8c/Project/yeast.nt /home/jill/bwa-0.5.8c/Project/yeast.fasta  > /home/jill/bwa-0.5.8c/Project/yeastbwaout.sai

#[bwa_aln] 17bp reads: max_diff = 2
#[bwa_aln] 38bp reads: max_diff = 3
#[bwa_aln] 64bp reads: max_diff = 4
#[bwa_aln] 93bp reads: max_diff = 5
#[bwa_aln] 124bp reads: max_diff = 6
#[bwa_aln] 157bp reads: max_diff = 7
#[bwa_aln] 190bp reads: max_diff = 8
#[bwa_aln] 225bp reads: max_diff = 9
#[bwa_aln_core] calculate SA coordinate... 0.15 sec
#[bwa_aln_core] write to the disk... 0.00 sec
#[bwa_aln_core] 1 sequences have been processed.

bwa samse /home/jill/bwa-0.5.8c/Project/yeast.nt  /home/jill/bwa-0.5.8c/Project/yeastbwaout.sai /home/jill/bwa-0.5.8c/Project/yeast.fasta  > /home/jill/bwa-0.5.8c/Project/yeastbwaout.sam 

#[bwa_aln_core] convert to sequence coordinate... 0.01 sec
#[bwa_aln_core] refine gapped alignments... 0.00 sec
#[bwa_aln_core] print alignments... 0.01 sec
#[bwa_aln_core] 1 sequences have been processed.

samtools faidx /home/jill/bwa-0.5.8c/Project/yeast.nt

samtools view -bt /home/jill/bwa-0.5.8c/Project/yeast.nt.fai /home/jill/bwa-0.5.8c/Project/yeastbwaout.sam > /home/jill/bwa-0.5.8c/Project/yeastbwaout.bam

#[samopen] SAM header is present: 17 sequences.

#index bam
samtools index /home/jill/bwa-0.5.8c/Project/yeastbwaout.bam

#Sorts bam
samtools sort /home/jill/bwa-0.5.8c/Project/yeastbwaout.bam /home/jill/bwa-0.5.8c/Project/yeastbwaoutsort 

samtools index /home/jill/bwa-0.5.8c/Project/yeastbwaoutsort.bam

samtools pileup -vcf /home/jill/bwa-0.5.8c/Project/yeast.nt.fai /home/jill/bwa-0.5.8c/Project/yeastbwaoutsort.bam | tee /home/jill/bwa-0.5.8c/Project/raw.txt | /home/jill/samtools-0.1.11/misc/samtools.pl varFilter -D100 > /home/jill/bwa-0.5.8c/Project/flt.txt
awk '($3=="*"&&$6>=50)||($3!="*"&&$6>=20)' /home/jill/bwa-0.5.8c/Project/flt.txt > /home/jill/bwa-0.5.8c/Project/final.txt
I made sure both my files have an equal line length, 80, formatted with fastx. Any insight would be very appreciated. Thanks so much!
mindlessbrain is offline   Reply With Quote
Old 12-12-2010, 12:39 PM   #2
drio
Senior Member
 
Location: 4117'49"N / 24'42"E

Join Date: Oct 2008
Posts: 323
Default

Try passing the actual fasta file instead of the of the index (fai) in the last command (samtools pileup).
__________________
-drd
drio is offline   Reply With Quote
Old 12-12-2010, 11:17 PM   #3
mindlessbrain
Junior Member
 
Location: qatar

Join Date: Dec 2010
Posts: 7
Default

Tried this:

Code:
samtools pileup -vcf /home/jill/bwa-0.5.8c/Project/yeast.nt /home/jill/bwa-0.5.8c/Project/yeastbwaoutsort.bam | tee /home/jill/bwa-0.5.8c/Project/raw.txt | /home/jill/samtools-0.1.11/misc/samtools.pl varFilter -D100 > /home/jill/bwa-0.5.8c/Project/flt.txt
awk '($3=="*"&&$6>=50)||($3!="*"&&$6>=20)' /home/jill/bwa-0.5.8c/Project/flt.txt > /home/jill/bwa-0.5.8c/Project/final.txt
Without the .fai . Still got empty pileup files.
mindlessbrain is offline   Reply With Quote
Old 12-13-2010, 03:40 AM   #4
mindlessbrain
Junior Member
 
Location: qatar

Join Date: Dec 2010
Posts: 7
Default

I've been trying to do the same thing with MAQ, that is get SNPS, and I'm getting an error there too. I'm not sure if the two are related or not though.

Code:
maq fasta2bfa /home/jill/maq/Project/yeast.nt /home/jill/maq/Project/yeast.nt.bfa
maq fasta2bfa /home/jill/maq/Project/yeast.fasta /home/jill/maq/Project/yeast.fasta.bfa
#-- 1 sequences have been converted.
maq match /home/jill/maq/Project/yeast.fasta.map /home/jill/maq/Project/yeast.nt.bfa /home/jill/maq/Project/yeast.fasta.bfa

#[ma_load_reads] loading reads...
#[ma_load_reads] set length of the first read as 4624380.
#[ma_load_reads] 1*2 reads loaded.
#[ma_longread2read] encoding reads... 2 sequences processed.
#[ma_match] set the minimum insert size as 4624381.
#[match_core] Total length of the reference: 12155026
#[match_core] round 1/3...
#[match_core] making index...
#[match_index_sorted] no reasonable reads are available. Exit!
ETA: I also tried the bam/sam script I have above with completely different files. Still, nothing.

Last edited by mindlessbrain; 12-13-2010 at 04:16 AM.
mindlessbrain is offline   Reply With Quote
Old 12-13-2010, 08:54 AM   #5
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

Why only one read processed?
swbarnes2 is offline   Reply With Quote
Old 12-13-2010, 10:04 PM   #6
mindlessbrain
Junior Member
 
Location: qatar

Join Date: Dec 2010
Posts: 7
Default

Quote:
Originally Posted by swbarnes2 View Post
Why only one read processed?
No idea. I was hoping someone could tell me. This is my first time using bwa/sam tools.
mindlessbrain is offline   Reply With Quote
Old 12-14-2010, 05:00 AM   #7
drio
Senior Member
 
Location: 4117'49"N / 24'42"E

Join Date: Oct 2008
Posts: 323
Default

Quote:
Originally Posted by mindlessbrain View Post
Tried this:

Code:
samtools pileup -vcf /home/jill/bwa-0.5.8c/Project/yeast.nt /home/jill/bwa-0.5.8c/Project/yeastbwaoutsort.bam | tee /home/jill/bwa-0.5.8c/Project/raw.txt | /home/jill/samtools-0.1.11/misc/samtools.pl varFilter -D100 > /home/jill/bwa-0.5.8c/Project/flt.txt
awk '($3=="*"&&$6>=50)||($3!="*"&&$6>=20)' /home/jill/bwa-0.5.8c/Project/flt.txt > /home/jill/bwa-0.5.8c/Project/final.txt
Without the .fai . Still got empty pileup files.
Try with /home/jill/maq/Project/yeast.fasta, instead of yeast.nt
__________________
-drd
drio is offline   Reply With Quote
Old 12-14-2010, 09:25 AM   #8
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

Well, this obviously won't work if you only have one read processed. Does your .sam file only have one non-header line in it? Something must be wrong with your fastq.
swbarnes2 is offline   Reply With Quote
Old 12-16-2010, 12:35 AM   #9
mindlessbrain
Junior Member
 
Location: qatar

Join Date: Dec 2010
Posts: 7
Default

I got it working! Somewhere along the line my fastq had been corrupted. So the commands I had above were fine, with the exception of swapping the reference for the .fai file.

Now I have an SNP file, from my last command:

Code:
samtools pileup -vcf /home/jill/bwa-0.5.8c/Project/yeast_mrna_genes.fasta /home/jill/bwa-0.5.8c/Project/yeastbwaoutsort.bam.bam | tee /home/jill/bwa-0.5.8c/Project/raw.txt | /home/jill/samtools-0.1.11/misc/samtools.pl varFilter -D100 > /home/jill/bwa-0.5.8c/Project/flt.txt
awk '($3=="*"&&$6>=50)||($3!="*"&&$6>=20)' /home/jill/bwa-0.5.8c/Project/flt.txt > /home/jill/bwa-0.5.8c/Project/final.txt
I've been trying to figure out what each column corresponds to. I know the first is the reference sequence name, then position of SNP, then SNPreference, SNPquery(?), after that though, no idea. The documentation has a slightly different format, except for the one they don't explain. =)

AB016599 1118 G R 215 215 36 81 a,,,a,,,.a,,a,,.,,,,,,,,,.aa,,,A.A.,.A....aA.AAAaAAa.,a,,a,,,,,,aaaaa,,,a,,a,,,.. CCACA0ACCBCC%CCCC%CC>CD%CCC@%CCCCCCCACCC?CCC%CACBBC:CCC%CCCCCCCBC<CCCC%CCC%C%CCAC
EF123147 124 T G 42 42 35 5 ^Fg^:g^Fg^Fg^Fg CC;@C
X84042 4 T A 36 36 25 3 ^:A^:A^:A ACC
X84042 5 C A 36 36 25 3 AAA ACC
mindlessbrain is offline   Reply With Quote
Old 12-16-2010, 05:29 AM   #10
drio
Senior Member
 
Location: 4117'49"N / 24'42"E

Join Date: Oct 2008
Posts: 323
Default

Explained here.
Also, since you are on it, try mpileup instead of pileup and compare results.
__________________
-drd
drio 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 09:05 AM.


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