Daniel:Notebook/GenomeMiner/2013-8-7: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Djacobse
No edit summary
>Djacobse
No edit summary
Line 10: Line 10:
  and deletions) needed to transform the read string into the reference string. Only present if SAM record is for
  and deletions) needed to transform the read string into the reference string. Only present if SAM record is for
   an aligned read.
   an aligned read.
Basically, these should be able to give me an idea of the number of mismatches.  I'm not really sure yet whether XM or NM is better, so I'm taking them both.  First I awked the data from the 3rd, and 12-20th columns.
'''Awk Command'''
awk '{print $3,\t,$12,\t,$13,\t,$14,\t,$15,\t,$16,\t,$17,\t,$18,\t,$19,\t,$20}' hcrp_onlymapped_hl155.sam >
/home/djacobse/hl155/hl155_3_12_andbeyond.txt
'''Perl Script'''
#!usr/bin/perl
use strict;
use warnings;
#should count the NM:i:digit column of the script and return
#the probe name, mismatches, and edit number
open READS, "hl155_3_12_andbeyond.txt";
while (<READS>) {
my($probe)=$_=~/(Probe_set\d_\d\d?\d?\d?\d?)/;
my($mismatch) = $_ =~/(XM:i:\d)/;
my($mismatch2)= $mismatch =~ /(\d)/ ;
my($editnum)= $_ =~ /(NM:i:\d)/;
my($editnum2) = $editnum =~ /(\d)/;
print $probe."\t".$mismatch2."\t".$editnum2."\n";
}

Revision as of 22:32, 7 August 2013

HL155

Back to Calendar

Continuing the analysis of the HL155 data. I looked through the information on the bowtie2 manual. The 12th column of a SAM output contains several possibly useful pieces of information:

From Bowtie2 manual

XM:i:<N> The number of mismatches in the alignment. Only present if SAM record is for an aligned read.
NM:i:<N> The edit distance; that is, the minimal number of one-nucleotide edits (substitutions, insertions
and deletions) needed to transform the read string into the reference string. Only present if SAM record is for
 an aligned read.

Basically, these should be able to give me an idea of the number of mismatches. I'm not really sure yet whether XM or NM is better, so I'm taking them both. First I awked the data from the 3rd, and 12-20th columns.

Awk Command

awk '{print $3,\t,$12,\t,$13,\t,$14,\t,$15,\t,$16,\t,$17,\t,$18,\t,$19,\t,$20}' hcrp_onlymapped_hl155.sam >
/home/djacobse/hl155/hl155_3_12_andbeyond.txt

Perl Script

#!usr/bin/perl
use strict;
use warnings;
#should count the NM:i:digit column of the script and return
#the probe name, mismatches, and edit number
open READS, "hl155_3_12_andbeyond.txt";
while (<READS>) {
my($probe)=$_=~/(Probe_set\d_\d\d?\d?\d?\d?)/;
my($mismatch) = $_ =~/(XM:i:\d)/;
my($mismatch2)= $mismatch =~ /(\d)/ ;
my($editnum)= $_ =~ /(NM:i:\d)/;
my($editnum2) = $editnum =~ /(\d)/;
print $probe."\t".$mismatch2."\t".$editnum2."\n";
}