SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
SOAP alignment format convert to SAM/BAM KevinLam Bioinformatics 31 01-10-2018 08:05 PM
SAM/BAM format to wiggle format pinki999 Bioinformatics 19 08-12-2015 12:35 AM
SAM to CUFFLINKS SAM format repinementer Bioinformatics 4 03-15-2012 08:53 AM
Looking process to convert gff3 format into ace format or sam format andylai Bioinformatics 1 05-17-2011 02:09 AM
anyone help me on bowtie format -> sam format! tninja Bioinformatics 2 04-25-2010 09:33 PM

Reply
 
Thread Tools
Old 08-05-2011, 01:04 AM   #281
Hena
Member
 
Location: Finland

Join Date: Nov 2009
Posts: 19
Default

Does anyone have nice perl snippet to calculate the end point of the read in a sam file?
Hena is offline   Reply With Quote
Old 10-20-2011, 04:26 AM   #282
pengchy
Senior Member
 
Location: China

Join Date: Feb 2009
Posts: 116
Default maq2sam-long problems

Dear all,

Quote:
Originally Posted by xmluo View Post
I used maq2sam-long to convert maq output to sam format, but all pairing information is missed in the results: MRNM is "*", and MPOS and ISIZE are 0. Can you recommend how to get these information? Thanks, Mei

I run into the same questions.

Another question is about the statistics for the results of out.map and out.bam.

1. how can I calculate the mapping rate?
95%(1903176 / 2000000) from maq or 95.68%(1820937 / 1903176) from the converted bam?

2. how can I calculate the unique mapped pair-end reads number?
I have conceived a method: according to the sam format spcification, the TAG "H0:i:1 H1:i:0", give the Number of perfect hits (H0) and Number of 1-di erence hits (H1), so we can deduce that one pair is uniquely mapped if and only if any one end's H0:i:1 or "H0:i:0 and H1:i:1", can it work?

3. how to get the MRNM, MPOS and ISIZE?
In my case, the ISIZE is normal. I found that the flag value in the second column is normal too. The mate read also given at another line, so we can fill the MRNM and MPOS columns accordingly. can it work?

Code:
##suppose we got the output file of "maq map" is out.map. then
##got out.sm
samtools-0.1.18/misc/maq2sam-long out.map  > out.sam
##got the sorted bam file out.bam
samtools-0.1.18/samtools view -ut reference.fa.fai out.sam | samtools  sort - out
##got the statistics results
##maq
maq-0.7.1/maq.pl statmap nohup.out > out.map.stat
samtools flagstat out.bam > out.bam.stat

-bash-3.2$ more out.map.stat 

-- == statmap report ==

-- # single end (SE) reads: 0
-- # mapped SE reads: 0 (/ 0 = NA%)
-- # paired end (PE) reads: 2000000
-- # mapped PE reads: 1903176 (/ 2000000 = 95.15%)
-- # reads that are mapped in pairs: 1756127 (/ 1903176 = 92.27%)
-- # Q>=30 reads that are moved to meet mate-pair requirement: 2962 (/ 1756127 = 0.16%)
-- # Q<30 reads that are moved to meet mate-pair requirement: 203102 (11.56%)

-bash-3.2$ more out.bam.stat 
1903176 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
1820937 + 0 mapped (95.68%:nan%)
1903176 + 0 paired in sequencing
951588 + 0 read1
951588 + 0 read2
1673888 + 0 properly paired (87.95%:nan%)
1738698 + 0 with itself and mate mapped
82239 + 0 singletons (4.32%:nan%)
1738698 + 0 with mate mapped to a different chr
1471758 + 0 with mate mapped to a different chr (mapQ>=5)
Any suggestions are appreciated, Thank you all.

pengchy
pengchy is offline   Reply With Quote
Old 10-21-2011, 05:48 PM   #283
pengchy
Senior Member
 
Location: China

Join Date: Feb 2009
Posts: 116
Default

Hi,
I have found the solution.
1. the mapping rate is 1820937/2000000

2. the unique mapped reads can be got from:
for the single end: "H0:i:1 H1:i:A" or "H0:i:0 H1:i:1" or "H0:i:0 H1:i:0" ##here, "A" stand for any number
for the pair end: any end has above tags.

3. Yes, it is.

Last edited by pengchy; 10-21-2011 at 05:51 PM. Reason: a mistake
pengchy is offline   Reply With Quote
Old 10-15-2012, 12:21 AM   #284
dtc
Junior Member
 
Location: ECNU

Join Date: Oct 2012
Posts: 3
Default

I have the same problem now, do you have any solution now?
Thanks Alexey
dtc is offline   Reply With Quote
Old 10-15-2012, 12:39 AM   #285
Hena
Member
 
Location: Finland

Join Date: Nov 2009
Posts: 19
Default

Quote:
Originally Posted by dtc View Post
I have the same problem now, do you have any solution now?
Thanks Alexey
I did my own version of this in the end using CIGAR string. However it does not take into account 'X', '=', 'N' or 'P' tags.

Code:
sub last_pos {
  my @read = @_;
  local $_;
  my ($match,$del) = (0,0);
  while ($read[$CIGAR] =~ m/(\d+)M/g) {
    $match += $1;
  }
  while ($read[$CIGAR] =~ m/(\d+)D/g) {
    $del += $1;
  }
  return $read[$POS] + $match + $del -1;
}
So counting the match/mismatch and deletion tags should give the number of nucleotides it will span in reference.
Hena is offline   Reply With Quote
Old 06-05-2013, 07:32 AM   #286
abi
Member
 
Location: Blacksburg

Join Date: May 2013
Posts: 18
Default SAM to BAM problem

Parse warning at line 27: mapped sequence without CIGAR
Parse error at line 27: sequence and quality are inconsistent
Aborted



I used BWASW for alignment. Any help ??
abi is offline   Reply With Quote
Old 06-05-2013, 07:49 AM   #287
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Quote:
Originally Posted by abi View Post
Parse warning at line 27: mapped sequence without CIGAR
Parse error at line 27: sequence and quality are inconsistent
Aborted



I used BWASW for alignment. Any help ??
Was the file truncated at that line? Otherwise it'd be helpful to see that line of the SAM file.
dpryan is offline   Reply With Quote
Old 06-05-2013, 08:01 AM   #288
abi
Member
 
Location: Blacksburg

Join Date: May 2013
Posts: 18
Default

No the file was not truncated at line 27. Here is the line:

0 chr4 96674301 60 227M * 0 0 TCTTTGAAACCAACGAGAACAAAGACACAACATACCAGAATCTCTGGGACACATTTAAAGCAGTGTGTAGAGGGAAATTTACAGCAGTAAATGCCCACAAGAGAATGCAGGAAAGATCTAAAATTGACACCCTAACATCAAAATTAAAAGAACTAGAGAAGCAAGAGCAATCACATTCAAAAGCTAGCAGAAGGCAAGAAATAACTAAGATCAGAGCAGAACTGAAG * NM:i:0 AS:i:227 XS:i:0


I am using bwa-0.7.4 for my alignment .

I must say that I used FASTA file for alignment and not FASTQ file. Here is a sample line from FASTA file:

> MySpecies 207
GGAGCCTTAAGGTTTGGTTATCGCCTTGTGTTTGTTTCTGGGGGTATCTGTGGGGTATGTGTTTCTGGCCATGTGTCTGTGTCTGTGTCTCTAGGCTGTCTTCTAGTCTCAGCTTGAGATCCACAGGCTTCAAGAGCTCAAGGGGGGAAAAGCCCAATTGTATATAAATTGTGAATGGGACTGATGCGTATGAGACAGGGAGGGTCT
abi is offline   Reply With Quote
Old 06-05-2013, 08:57 AM   #289
sbaheti
Member
 
Location: Rochester

Join Date: Jul 2010
Posts: 12
Default

I have seen similar issues with BWA 0.7 greater, same command runs absolutely fine with 0.7 version of BWA.
sbaheti is offline   Reply With Quote
Old 06-05-2013, 09:38 AM   #290
narain
Member
 
Location: Washington DC

Join Date: Aug 2011
Posts: 78
Thumbs up

@abi , lh3, sbaheti, dpryan, Hena and rest

I have the answer to you. There is a bug in newer versions of BWA. It does not produces right SAM files for those sequences which have their Fasta/Fastq files starting with > and a space. As an example:

This will not produce right SAM files
> My sequence

This will produce right SAM files
>My sequence


Notice the difference in space after '>'

Hope this helps. You can read more about me at www.tinyurl.com/abinarain

Cheers
Narain

Last edited by narain; 06-05-2013 at 09:40 AM.
narain is offline   Reply With Quote
Old 01-21-2014, 10:37 AM   #291
rpauly
Member
 
Location: Atlanta

Join Date: Apr 2011
Posts: 32
Default

[samopen] SAM header is present: 84 sequences.
[sam_read1] reference '121337031' is recognized as '*'.
Parse warning at line 7932555: mapped sequence without CIGAR
Parse error at line 7932555: sequence and quality are inconsistent
Aborted (core dumped)


I am getting this error... tried the suggestion on changing the MDI in the blast2sam.pl but it is still not working!

Any suggestions will be appreciated.
rpauly is offline   Reply With Quote
Old 01-21-2014, 10:42 AM   #292
rpauly
Member
 
Location: Atlanta

Join Date: Apr 2011
Posts: 32
Default

[samopen] SAM header is present: 84 sequences.
[sam_read1] reference '121337031' is recognized as '*'.
Parse warning at line 7932555: mapped sequence without CIGAR
Parse error at line 7932555: sequence and quality are inconsistent
Aborted (core dumped)

The suggestion to change MDI in blast2sam pl did not work! Any suggestions will be appreciated.
rpauly is offline   Reply With Quote
Old 01-22-2014, 12:25 AM   #293
TiborNagy
Senior Member
 
Location: Budapest

Join Date: Mar 2010
Posts: 329
Default

Blast2Sam won't work, because it is an outdated script (I have tried it several times without success). Try to align your sequences with bwa/bowtie.
TiborNagy 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 01:29 AM.


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