SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Perl script to convert BAM to BED? rebrendi Bioinformatics 3 03-22-2012 10:34 PM
just perl script semna Bioinformatics 3 07-02-2011 08:42 AM
vcftools perl script weiyulin Bioinformatics 6 12-09-2010 02:13 PM
perl script bioenvisage Bioinformatics 0 02-01-2010 07:23 AM
Perl script bioenvisage Bioinformatics 4 01-28-2010 12:25 PM

Reply
 
Thread Tools
Old 01-26-2010, 11:55 AM   #1
bioenvisage
Member
 
Location: it

Join Date: Oct 2009
Posts: 40
Default perl script

Hai .I need a perl script can anybody help me ...I did an alignment of illumina paired reads,I need a script to extract only those reads in which both mates are not placed.
example of the file which i need to extract

HWI-EAS373:1:1:10:1119#0/1 GGAGCATTGCGACATTTAGCTNATTCAGCGGCATGGAAGG abbbbbbbbbaaa^aaa`aa^D[_`a_[_^__Z_a\G[__ NF
HWI-EAS373:1:1:10:1119#0/2 GGAATAAGCAGAGTGAGCATAAGAGAAGATGNCTTCATGC a^WW_^`a`W^YXS[T\^M[YH[\_VY_]]^D[^^YPY]` NF
HWI-EAS373:1:1:10:948#0/1 CTGGTCTTCCTCGATTCTCCCNATGGAGCCACTTTCTTCA abb`bbabbbbaa[a_abba^DP__[TX_^UIT[[___^U NF
HWI-EAS373:1:1:11:1578#0/1 GGCAGGGTGAGGCTCCAGTCANTTGTCACACCCTTAACGG aaXababUbabb`\``[\^]VDWXOV]__]_`\MSKKa\Y NF

The last two reads are NF (no match found), but their mate pair is aligned , so i need to extract only the reads in which both mates are not placed (NF)..

thanks
bioenvisage is offline   Reply With Quote
Old 01-26-2010, 01:21 PM   #2
jnfass
Member
 
Location: Davis, CA

Join Date: Aug 2008
Posts: 88
Default perl one-liners to pull pairs of unmapped reads

I often use perl on the *nix command line to solve problems like this. For this example (assuming the file you gave an example of is named "export.txt", you could:

cat export.txt | perl -ane 'print if ($F[3] eq "NF")' | cut -f1 -d\# | sort | uniq -c | grep " 2 " | perl -ne 'm/(HWI.*)$/; print ">".$1."\n"' > pairedIds.txt

That should give you the id's (without the "#0/1" or "/2", but with prepended ">") of the reads that appear twice as "NF" ... you can then pull these from the two s_?_?_sequence.txt files with:

grep -A3 -F -f pairedIds.txt s_?_1_sequence.txt | grep -v ^-- > F.txt
grep -A3 -F -f pairedIds.txt s_?_2_sequence.txt | grep -v ^-- > R.txt

These are big grep's; they'll take a while. Note that grep can introduce separator lines of only two dashes when it finds consecutive blocks matching the pattern; the grep -v is intended to take those.

Also, the above is 'off the cuff' ... you might have to debug a little
jnfass is offline   Reply With Quote
Old 01-27-2010, 12:08 AM   #3
simonandrews
Simon Andrews
 
Location: Babraham Inst, Cambridge, UK

Join Date: May 2009
Posts: 871
Default

Assuming you want the full lines for the entries which were not found, and that the sequence pairs come immediately after one another then something like this should work:

Code:
#!/usr/bin/perl
use warnings;
use strict;

my $last_line;
my $last_id = '';

while(<>) {

  # Ignore any lines for sequences which were found
  next unless (/NF$/);

  # Extract the id minus the strand
  my $id;
  if (/^([^\/]+)/) {
    $id = $1;
  }
  else {
    die "Couldn't get id from $_";
  }

  # If this line has the same id as the last one then
  # both were not found so print them out
  if ($id eq $last_id) {
    print $last_line,$_;
  }

  # If not then store this line in case the next one
  # forms a pair which we'll print later.
  else {
    $last_line = $_;
    $last_id = $id;
  }
}
This will read any data passed on stdin or filenames passed as an argument and will output on stdout.
simonandrews is offline   Reply With Quote
Old 01-27-2010, 01:07 AM   #4
bioenvisage
Member
 
Location: it

Join Date: Oct 2009
Posts: 40
Default

hi thanks for the scripts ..the second script from simon andrews is very fast ..
bioenvisage is offline   Reply With Quote
Old 02-01-2010, 05:10 AM   #5
bioenvisage
Member
 
Location: it

Join Date: Oct 2009
Posts: 40
Default perl script

................

Last edited by bioenvisage; 02-01-2010 at 11:53 PM.
bioenvisage is offline   Reply With Quote
Old 02-01-2010, 08:11 AM   #6
simonandrews
Simon Andrews
 
Location: Babraham Inst, Cambridge, UK

Join Date: May 2009
Posts: 871
Default

I see you've posted another 'perl script' thread on the forum with another related question. I don't mean to sound harsh, but this isn't a place to come to get work done for you - though I'm sure there are people here who have consultancy services.

I don't mind proving occasional scripts when people get stuck - and I'm more than happy to help sort out problems with scripts people have tried writing for themselves, but there's a limit to how much you can expect without at least showing what you've tried to sort this out yourself.

If you're going to be doing this sort of work then I'd strongly recommend finding a local Perl course. I'm sure you'd find that a lot more productive than repeatedly posting to get canned answers. Using scripts you got from a forum without completely understanding what they do isn't going to be able to give you much confidence in your analysis results.

If you've tried to modify the previously posted script and couldn't get it to work please feel free to post that and I'm sure you'll get help to make it work.
simonandrews 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 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO