Mercurial > repos > tyty > structurefold
changeset 105:abe41dc167f1 draft
Uploaded
author | tyty |
---|---|
date | Thu, 19 Mar 2015 17:43:24 -0400 |
parents | 60278e5df493 |
children | c1c1d38ec295 |
files | ._get_reads get_reads/._.DS_Store get_reads/._get_read.py get_reads/._get_read.xml get_reads/._read_file.py get_reads/get_read.py get_reads/get_read.xml get_reads/read_file.py get_reads/read_file.pyc |
diffstat | 9 files changed, 147 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/get_reads/get_read.py Thu Mar 19 17:43:24 2015 -0400 @@ -0,0 +1,80 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +import sys +from Bio import SeqIO +import os +from read_file import * +import random +import string + +fasta_file = sys.argv[1] +map_file = sys.argv[2] +result_file = sys.argv[3] + +syspathrs = os.getcwd() + +os.system("samtools view -F 0xfff "+map_file+"|cut -f 3,4 > "+syspathrs+"map_info.txt") + +fasta_sequences = SeqIO.parse(open(fasta_file),'fasta'); +length_seq = {}; +for seq in fasta_sequences: + nuc = seq.id; + length_seq[nuc] = len(seq.seq.tostring()); + + + +mapping = {} +transcripts = [] + +f = open(syspathrs+"map_info.txt"); +for aline in f.readlines(): + tline = aline.strip(); + tl = tline.split('\t'); + if tl[0].strip() not in transcripts: + transcripts.append(tl[0].strip()); + mapping[tl[0].strip()] = []; + + mapping[tl[0].strip()].append(tl[1].strip()); + +distribution = {}; +coverage = {}; +for transcript in length_seq: + distribution[transcript] = []; + for i in range(0, length_seq[transcript]): + distribution[transcript].append(0); + sum_count = float(0); + if transcript in mapping: + for j in range(0, len(mapping[transcript])): + index = mapping[transcript][j]; + #count = reads[mapping[transcript][j][0]]; + sum_count = sum_count + 1; + distribution[transcript][int(index)-1] = distribution[transcript][int(index)-1] + 1; + coverage[transcript] = float(sum_count)/float(length_seq[transcript]); + else: + coverage[transcript] = 0 + + + + + +h = file(result_file, 'w') +for transcript in length_seq: + h.write(transcript); + h.write('\n') + for i in range(0, length_seq[transcript]): + h.write(str(distribution[transcript][i])) + h.write('\t') + h.write('\n') + h.write('\n') + +#os.system("rm -r "+syspathrs) + + + +f.close(); +h.close() + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/get_reads/get_read.xml Thu Mar 19 17:43:24 2015 -0400 @@ -0,0 +1,46 @@ +<tool id="get_read_pipeline" name="Get RT Stop Counts" version="1.0"> + <description>derives the reverse transcriptase (RT) stop count on each nucleotide from a mapped file provided by the Iterative Mapping module</description> + <command interpreter="python">get_read.py $lib_file $map_file $output </command> + <requirements> + <requirement type="package" version="1.61">biopython</requirement> + <requirement type="package" version="1.7.1">numpy</requirement> + <requirement type="package" version="0.1.18">samtools</requirement> + </requirements> + <inputs> + <param name="lib_file" type="data" format="fasta" label="Reference genome/transcriptome"/> + <param name="map_file" type="data" format="bam" label="Mapped file"/> + </inputs> + <outputs> + <data name="output" format="txt"/> + </outputs> + <tests> + <test> + <param name="lib_file" value="test.bam" /> + <param name="map_file" value="com_rna.txt" /> + <output name="output" file="get_RT_stop_test.out" /> + </test> + </tests> + <help> + + +**Function** + +Get RT Stop Counts derives the RT stop counts on each nucleotide of each transcript in the reference transcriptome based on a mapped file (.bam), typically the output from the Iterative Mapping module. + +----- + +**Input**: + +* 1. A mapped (.bam) file from Bowtie (or any other mapping program) +* 2. Reference library sequences (fasta) used to map the reads to + +----- + +**Output**: + +A text file with reverse transcription stop counts mapped to each nucleotide (RTSC file) + + + + </help> +</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/get_reads/read_file.py Thu Mar 19 17:43:24 2015 -0400 @@ -0,0 +1,21 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +import sys + + + +def read_t_file(in_file): + f = open(in_file); + result = []; + for aline in f.readlines(): + temp = []; + tline = aline.strip(); + tl = tline.split('\t'); + for i in range(0, len(tl)): + temp.append(tl[i].strip()); + result.append(temp); + f.close(); + return result; + +