-
Notifications
You must be signed in to change notification settings - Fork 0
/
Snakefile
79 lines (71 loc) · 2.59 KB
/
Snakefile
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
import os
from script.snakefile_functions import verify_output_directory
configfile: "config.yaml"
config["output_directory"] = verify_output_directory(config["output_directory"])
proteome_list = os.listdir(config["proteomes"])
blast_to_use = "blastx" if config["is_transcriptome"] == "yes" else "blastp"
#
rule all:
input:
config["output_directory"] + "Statistics.txt"
# Take all proteome in given proteome folder path and create
# the corresponding database with makeblastdb.
rule makeblastdb:
input:
config["proteomes"] + "/{proteome_list}"
output:
config["output_directory"] + "blastdb/{proteome_list}"
shell:
"""
ln -s {input} {output}
makeblastdb -in {output} -dbtype prot
"""
rule blastx:
input:
proteomes = config["output_directory"] + "blastdb/{proteome_list}",
transcriptome = config["transcriptome"]
output:
config["output_directory"] + "blast_result/{proteome_list}_blastx.txt"
threads:
config["threads"]
params:
cutoff = config["evalue"],
blast_to_use = blast_to_use
shell:
"""
{params.blast_to_use} -query {input.transcriptome} -db {input.proteomes} \\
-num_threads {threads} -max_hsps 1 -outfmt \"6 std stitle\" \\
-evalue {params.cutoff} > {output}
"""
# generate a file listing all transcripts.
rule transcript_list:
input:
transcriptome = config["transcriptome"]
output:
config["output_directory"] + "transcript_list.txt"
shell:
"grep \"^>\" {input} | cut -d ' ' -f1 | tr -d \">\" | sort | uniq > {output}"
rule transcript_to_annotation:
input:
trans_list=config["output_directory"] + "transcript_list.txt",
blast=config["output_directory"] + "blast_result/{proteome_list}_blastx.txt"
output:
config["output_directory"] + "{proteome_list}_transcript_to_annotation.txt"
shell:
# bash script/trans_to_ann.sh {input.trans_list} {input.blast} {output}
"""
cat {input.blast} | cut -f1,13 | sort -u -k1,1 > {output}
cat {output} | cut -f1 > {output}tmp.txt
res=$( cat {output}tmp.txt {input.trans_list} | sort | uniq -u )
for i in $res; do
echo -e "${{i}}\tnone" >> {output}
done
rm {output}tmp.txt
"""
rule compute_statistics:
input:
expand(config["output_directory"] + "{proteome_list}_transcript_to_annotation.txt", proteome_list = proteome_list)
output:
config["output_directory"] + "Statistics.txt"
shell:
"bash script/compute_statistics.sh {output} {input}"