Daniel:Notebook/GenomeMiner/2013-7-31

From ZhangLabWiki
Jump to navigation Jump to search

HL155[edit]

Back to Calendar

Started with the sam file of only mapped reads. Then printed the third column into text file.

Awk third column from SAM

awk '{print $3}' hrcp_onlymapped.sam > test.txt

Next I used perl to count the reads using this script:

#!usr/bin/perl
#reads column 3 entries from hrcp_onlymapped_hl155.sam
#finds probe entries "Probe_setn_m" and adds to probe counts of variable $probearray
my @probearray = ([0]x19518,[0]x20000,[0]x18850,[0]x19301);
my $set = 0;
open READS, "test.txt";
while (<READS>) {
if (/Probe/) {
$set = 1 if (/set1/);
$set = 2 if (/set2/);
$set = 3 if (/set3/);
$set = 4 if (/set4/);
my ($probe,$junk,$num)=split /_/, $_;
chomp($num);
++$probearray[$set][$num];
}
}
for (1..19517) {
print "$probearray[1][$_],";
}
print "$probearray[1][19518]\n";
for (1..19999) {
print "$probearray[2][$_],";
}
print "$probearray[2][20000]\n";
for (1..18849) {
print "$probearray[3][$_],";
}
print "$probearray[3][18850]\n";
for (1..19300) {
print "$probearray[4][$_],";
}
print "$probearray[4][19301]\n";

Copied probecount file back to my computer for processing in MATLAB.

SCP to laptop

%% scp <remote user>@<remote machine>:<remote path> <local file>
scp djacobse@132.239.135.41:/home/djacobse/HL155/130628_HL155/probecount.csv probecount.csv