Mercurial > repos > artbio > small_rna_maps

comparison test-data/test/small_rna_map.bash @ 0:6d48150495e3 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_maps commit d4d8106d66b65679a1a685ab94bfcf99cdb7b959
author artbio
date Mon, 24 Jul 2017 06:28:45 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:6d48150495e3
1 #! /bin/bash
2 # Execution command : bash small_rna_map.bash <sam file>
3 samtools view -h -o input.sam ../test-data/input.bam
4 samfile="input.sam"
5 cut -f2,3,4,5,6 ../test-data/output.tab > test-data.dat
6 grep -e "^@SQ" $samfile > header.tab
7 cat $samfile | perl -ne 'print unless /^@/ or ( (split /\t/)[1] != 0 and (split /\t/)[1] != 16 )' > output1.tab
8
9 cat << EOF > pyscript.py
10 lendic = {}
11 for line in open("header.tab"):
12 fields = line[:-1].split()
13 lendic[fields[1][3:]] = fields[2][3:]
14 readdic = {}
15 for line in (open("output1.tab")):
16 fields = line[:-1].split()
17 key = "-".join([fields[2],fields[3],"F" if fields[1]=="0" else "R"])
18 try:
19 readdic[key] += 1
20 except KeyError:
21 readdic[key] = 1
22 Out = open("output.tab", "w")
23 Out.write("Chromosome\tChrom_length\tCoordinate\tNbr_reads\tPolarity\n")
24 for key in sorted(readdic):
25 fields = key.split("-")
26 Chromosome, Coordinate, Polarity = fields[0], fields[1], fields[2]
27 Chrom_length = lendic[Chromosome]
28 Nbr_reads = str(readdic[key])
29 Out.write("\t".join([Chromosome, Chrom_length, Coordinate, Nbr_reads, Polarity]) )
30 Out.write("\n")
31 Out.close
32 EOF
33
34 python ./pyscript.py
35 rm output1.tab header.tab pyscript.py
36 sort -k 1,1 -k 3,3n -k 5,5 output.tab -o output.tab
37 if ! diff -q output.tab test-data.dat &>/dev/null; then
38 >&2 echo "different"
39 else
40 >&2 echo "Test passed"
41 fi
42
43 rm output.tab input.sam test-data.dat