View Single Post
Old 05-24-2012, 03:47 AM   #21
cedance
Senior Member
 
Location: Germany

Join Date: Feb 2011
Posts: 108
Default

I map it on to the genome. I do not use picard. I wrote a perl script using Bio:B::Sam, with which its very easy to obtain all pairs and their inner distance. I take only uniquely mapped pairs and where quality is >= 20 to compute the -r and --mate-std-dev parameter. It goes something like this:

Quote:
# $opt_i is your input bam file, bam.bai must also exist in the same name under the same directory
my $sam = Bio:B::Sam->new( -bam => $opt_i );
# get all chromosome ids
my @targets = $sam->seq_ids;
# for each chromosome id
foreach my $seqid (@targets ) {
my $segment = $sam->segment( -seq_id => $seqid ); # get all reads that match
# get iterator to loop over all pairs
my $iterator = $segment->features(-type => 'read_pair', -iterator => 1 );
while ( my $pair = $iterator->next_seq ) {
# fetch 1 pair at a time
my ( $first_mate, $second_mate ) = $pair->get_SeqFeatures;
# check conditions to skip to next pair or not
next if( !defined( $second_mate) );
next if( $first_mate->get_tag_values( "XT" ) ne "U" );
next if( $first_mate->qual < 20 );
# conditions cleared? get inner distance
my $idist = ($second_mate->start - 1) - ($first_mate->end + 1) + 1;

# here, you can just push it to an array and after you get out of the for-loop compute the quantiles and the IQ and the mean and SD.
....
....
}
}

# compute mean and SD by computing IQ and filtering those inner distances that are within Q1 - (Q3-Q1)*2 to Q3 + (Q3-Q1)*2, where Q3-Q1 is the IQ (inter-quartile range).
cedance is offline   Reply With Quote