perl script

[bioenvisage]

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

[jnfass]

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

[simonandrews]

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:

#!/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.

[bioenvisage]

hi thanks for the scripts ..the second script from simon andrews is very fast ..

[bioenvisage]

perl script

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

[simonandrews]

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.