[TiborNagy]

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.

[rpauly]

[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]

[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.

[narain]

@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

[sbaheti]

I have seen similar issues with BWA 0.7 greater, same command runs absolutely fine with 0.7 version of BWA.

[abi]

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

[dpryan]

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.

[abi]

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 ??

[Hena]

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.

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.

[dtc]

I have the same problem now, do you have any solution now?
Thanks Alexey

[pengchy]

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.

[pengchy]

maq2sam-long problems

Dear all,

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-dierence 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?

##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

[Hena]

Does anyone have nice perl snippet to calculate the end point of the read in a sam file?

[Bruins]

Dear all,

I'm had the same problem as oudacontrol above and Jeckow earlier in this thread.

Trying to get samtools view to extract one chromosome from a sorted and indexed bam file fails:

[sam_header_read2] 84 sequences loaded.
[main_samview] random alignment retrieval only works for indexed BAM files.

It took me a long while to figure out what caused it, so I thought I'd post my simple solution for poor sobs like me...

samtools sort in.bam in.sorted

samtools index in.sorted.bam

samtools view -bh -t human_g1k_v37.fa.fai -o in.sorted.chr9.bam in.sorted.bam 9

As it turned out, I should not have used the -t option. Without it, view finished as it should have, leaving me with a nice chr9-only bam file! (at least, it's a 279M file, I'll have to check later if trackster or igv will load it )

Cheers,
Bruins

[jnfass]

@oudacontrol: You should post the command you used. It seems like the command may have been mis-formatted, but there's no way to tell if you don't include it.