copangraph/copangraph.smk at main · korem-lab/copangraph · GitHub

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
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
from os.path import join
import pandas as pd
import sys


config_path = sys.argv[sys.argv.index('--configfile')+1]
sample_dir = config['sample_dir']
out_dir = config['out_dir']
assembler = config['assembler']
metadata = config['metadata']

# get samples
metadata_csv =  pd.read_csv(join(sample_dir, metadata))
SAMPLES = [str(s) for s in metadata_csv['SAMPLES']]
print('INPUT SAMPLES:')
for e in SAMPLES:
    print(e)

rule all:
    input:
        join(out_dir, 'copangraph', config['graph_name'], config['graph_name'] + '.gfa')

rule extension:
    input:
        expand(join(out_dir, 'extended_contigs/{sample}.pe_ext.fasta'), sample=SAMPLES)

rule assemble_samples_megahit:
    input:
        r1=join(sample_dir, '{sample}_1.fastq.gz'),
        r2=join(sample_dir, '{sample}_2.fastq.gz')
    output:
        fl=join(out_dir, assembler, '{sample}/final.contigs.fa'),
        dir=directory(join(out_dir, assembler, '{sample}'))
    log:
        join(out_dir, assembler, '{sample}/assemble_samples_megahit_log.txt')
    threads:
        8
    shell:
        'megahit -t {threads} -1 {input.r1} -2 {input.r2} -f -o {output.dir} > {log} 2>&1'

rule assemble_samples_metaspades:
    input:
        r1=join(sample_dir, '{sample}_1.fastq.gz'),
        r2=join(sample_dir, '{sample}_2.fastq.gz')
    output:
        fl=join(out_dir, assembler, '{sample}/contigs.fasta'),
        dir=directory(join(out_dir, assembler, '{sample}'))
    log:
        join(out_dir, assembler, '{sample}/assemble_samples_metaspades_log.txt')
    threads:
        8
    shell:
        'metaspades.py -t {threads} -1 {input.r1} -2 {input.r2} -f -o {output.dir} > {log} 2>&1'

rule map_reads:
    input:
        contigs=join(out_dir, assembler,'{sample}/final.contigs.fa' if assembler == 'megahit' else '{sample}/contigs.fasta'),
        r1=join(sample_dir, '{sample}_1.fastq.gz'),
        r2=join(sample_dir, '{sample}_2.fastq.gz')
    output:
        join(out_dir, 'read_mappings/{sample}/mapping.bam')
    log:
        build=join(out_dir, 'read_mappings/{sample}/map_reads_build_log.txt'),
        map=join(out_dir, 'read_mappings/{sample}/map_reads_map_log.txt')
    threads:
        8
    shell:
        'bowtie2-build --threads {threads} {input.contigs} {out_dir}/read_mappings/{wildcards.sample}/contigs > {log.build}  2>&1 && '
        'bowtie2 --threads {threads} -x {out_dir}/read_mappings/{wildcards.sample}/contigs -1 {input.r1} -2 {input.r2} 2> {log.map} | samtools view -@ {threads} -bS -h - > {out_dir}/read_mappings/{wildcards.sample}/mapping.bam'

rule sort_bam:
    input:
        join(out_dir, 'read_mappings/{sample}/mapping.bam')
    output:
        join(out_dir, 'read_mappings/{sample}/sorted_mapping.bam')
    threads:
        8
    log:
        join(out_dir, 'read_mappings/{sample}/sort_bam_log.txt')
    shell:
        'samtools sort -n -@ {threads} -o {output} {input} '

rule run_extension:
    input:
        bam=join(out_dir, 'read_mappings/{sample}/sorted_mapping.bam'),
        contigs=join(out_dir, assembler, '{sample}/final.contigs.fa' if assembler == 'megahit' else '{sample}/contigs.fasta')
    params:
        pe_only='--pe-only' if config.get('EXTENSION_pe_only', False) else '',
        depth=str(config.get('EXTENSION_depth', 10))
    threads:
        8
    output:
        join(out_dir, 'extended_contigs/{sample}.pe_ext.fasta')
    log:
        join(out_dir, 'extended_contigs/{sample}.run_extension_log.txt')
    shell:
        'bin/release/extension -i {input.contigs} -b {input.bam} -m {params.depth} {params.pe_only} -o ' + join(out_dir, 'extended_contigs') + ' -n {wildcards.sample} > {log} 2>&1 '


rule build_copangraph_config:
    input:
        contigs=expand(join(out_dir, 'extended_contigs/{sample}.pe_ext.fasta'), sample=SAMPLES)
    output:
        join(out_dir, 'copangraph', config['graph_name'], 'config.ini'),
        join(out_dir, 'copangraph', config['graph_name'], 'sample_list.txt')
    params:
        config_path
    script:
        'scripts/build_copangraph_config.py'


rule run_copangraph:
    input:
        join(out_dir, 'copangraph', config['graph_name'], 'config.ini'),
        join(out_dir, 'copangraph', config['graph_name'], 'sample_list.txt')
    output:
        join(out_dir, 'copangraph', config['graph_name'], config['graph_name'] + '.gfa')
    threads:
        8
    log:
        join(out_dir, 'copangraph', config['graph_name'], 'run_copangraph_log.txt')
    shell:
        'bin/release/copangraph -i {input} > {log} 2>&1'