SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Make my own blast DB detq182 Bioinformatics 9 04-18-2012 01:45 PM
Zygosity difference oyvindbusk Bioinformatics 3 01-31-2012 05:56 AM
Can't make Maq nightsun Bioinformatics 1 12-30-2011 08:24 PM
Make Clean and Make all not working qnc Bioinformatics 27 10-21-2011 10:17 AM
Can't make Maq RockChalkJayhawk Bioinformatics 4 08-24-2009 10:32 AM

Reply
 
Thread Tools
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
Old 05-24-2012, 04:32 AM   #22
EGrassi
Member
 
Location: Turin, Italy

Join Date: Oct 2010
Posts: 66
Default

I see, thanks. I did this evaluation with picard (and separately with awk on the bam file, should be similar to your perl) but with alignments done on the transcriptome because it seemed the sensible thing to do to me...maybe I was utterly wrong...ASAP I will perform some tests and report back here.

(the $idist formula could be cleaned removing some -1/+1 I believe , as in $second_mate->start - $first_mate->end -1)
EGrassi is offline   Reply With Quote
Old 05-24-2012, 04:38 AM   #23
cedance
Senior Member
 
Location: Germany

Join Date: Feb 2011
Posts: 108
Default

EGrassi, About the cleanup, Yes! However, its just to remember the next time I look at the code that I am computing the coordinates of the junctions at open intervals .

I just dint want to introduce any variability by mapping and determining inner distance and standard deviations from transcriptome. I don't think picard estimates inner distance this way. I remember using picard initially and I was not satisfied with the mapping.
cedance is offline   Reply With Quote
Old 05-25-2012, 04:43 AM   #24
glados
Member
 
Location: Aperture Science

Join Date: Mar 2012
Posts: 59
Default

Yes it's a good idea to calculate your mean inner distance and I would not trust the flag "properly paired" tophat gives. See this thread.
glados 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 11:39 PM.


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