Newer
Older
"""General purpose RNA-Seq analysis pipeline developed by the Zavolan Lab"""
BIOPZ-Gypas Foivos
committed
import os
BIOPZ-Gypas Foivos
committed
# Get sample table
samples_table = pd.read_csv(
config['samples'],
header=0,
index_col=0,
comment='#',
engine='python',
sep="\t",
)
BIOPZ-Katsantoni Maria
committed
def get_sample(column_id, search_id=None, search_value=None):
if search_id:
if search_id == 'index':
return str(samples_table[column_id][samples_table.index == search_value][0])
else:
return str(samples_table[column_id][samples_table[search_id] == search_value][0])
else:
return str(samples_table[column_id][0])
localrules: start, finish, rename_star_rpm_for_alfa, prepare_multiqc_config
BIOPZ-Katsantoni Maria
committed
if cluster_config:
os.makedirs(
os.path.join(
os.getcwd(),
os.path.dirname(cluster_config['__default__']['out']),
),
BIOPZ-Katsantoni Maria
committed
exist_ok=True)
BIOPZ-Gypas Foivos
committed
# Include subworkflows
include: os.path.join("workflow", "rules", "paired_end.snakefile.smk")
include: os.path.join("workflow", "rules", "single_end.snakefile.smk")
BIOPZ-Katsantoni Maria
committed
"""
Rule for collecting outputs
"""
multiqc_report = os.path.join(
config['output_dir'],
"multiqc_summary"),
bigWig = expand(
config["output_dir"],
"samples",
"{sample}",
"bigWig",
"{unique_type}",
"{sample}_{unique_type}_{strand}.bw"),
BIOPZ-Katsantoni Maria
committed
sample=pd.unique(samples_table.index.values),
strand=["plus", "minus"],
unique_type=["Unique", "UniqueMultiple"]),
salmon_merge_genes = expand(
os.path.join(
config["output_dir"],
"summary_salmon",
"quantmerge",
"genes_{salmon_merge_on}.tsv"),
salmon_merge_on=["tpm", "numreads"]),
salmon_merge_transcripts = expand(
os.path.join(
config["output_dir"],
"summary_salmon",
"quantmerge",
"transcripts_{salmon_merge_on}.tsv"),
salmon_merge_on=["tpm", "numreads"]),
rule start:
'''
Get samples
'''
input:
reads = lambda wildcards:
BIOPZ-Katsantoni Maria
committed
expand(
pd.Series(
samples_table.loc[wildcards.sample, wildcards.mate]
).values)
output:
reads = os.path.join(
config["output_dir"],
"samples",
"{sample}",
"start",
"{sample}.{mate}.fastq.gz")
log:
stderr = os.path.join(
config["log_dir"],
"samples",
"{sample}",
"start_{sample}.{mate}.stderr.log"),
stdout = os.path.join(
config["log_dir"],
"samples",
"{sample}",
"start_{sample}.{mate}.stdout.log")
BIOPZ-Katsantoni Maria
committed
"(cat {input.reads} > {output.reads}) \
1> {log.stdout} 2> {log.stderr} "
rule fastqc:
'''
A quality control tool for high throughput sequence data
'''
input:
reads = os.path.join(
config["output_dir"],
"samples",
"{sample}",
"start",
"{sample}.{mate}.fastq.gz")
output:
outdir = directory(
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
config["output_dir"],
"samples",
"{sample}",
"fastqc",
"{mate}"))
threads: 2
singularity:
"docker://zavolab/fastqc:0.11.9-slim"
log:
stderr = os.path.join(
config["log_dir"],
"samples",
"{sample}",
"fastqc_{mate}.stderr.log"),
stdout = os.path.join(
config["log_dir"],
"samples",
"{sample}",
"fastqc_{mate}.stdout.log")
shell:
"(mkdir -p {output.outdir}; \
fastqc --outdir {output.outdir} {input.reads}) \
1> {log.stdout} 2> {log.stderr}"
BIOPZ-Iborra de Toledo Paula
committed
BIOPZ-Gypas Foivos
committed
rule create_index_star:
BIOPZ-Katsantoni Maria
committed
"""
Create index for STAR alignments
"""
input:
genome = lambda wildcards:
BIOPZ-Katsantoni Maria
committed
get_sample(
'genome',
search_id='organism',
search_value=wildcards.organism),
BIOPZ-Katsantoni Maria
committed
get_sample(
'gtf',
search_id='organism',
search_value=wildcards.organism)
BIOPZ-Katsantoni Maria
committed
output:
chromosome_info = os.path.join(
config['star_indexes'],
"{organism}",
"{index_size}",
"STAR_index",
BIOPZ-Katsantoni Maria
committed
"chrNameLength.txt"),
chromosomes_names = os.path.join(
config['star_indexes'],
"{organism}",
"{index_size}",
"STAR_index",
BIOPZ-Katsantoni Maria
committed
"chrName.txt")
params:
output_dir = os.path.join(
config['star_indexes'],
"{organism}",
"{index_size}",
BIOPZ-Katsantoni Maria
committed
"STAR_index"),
outFileNamePrefix = os.path.join(
config['star_indexes'],
"{organism}",
BIOPZ-Katsantoni Maria
committed
"STAR_index/STAR_"),
sjdbOverhang = "{index_size}"
BIOPZ-Katsantoni Maria
committed
BIOPZ-Katsantoni Maria
committed
BIOPZ-Katsantoni Maria
committed
stderr = os.path.join(
BIOPZ-Katsantoni Maria
committed
"{organism}_{index_size}_create_index_star.stderr.log"),
stdout = os.path.join(
config['log_dir'],
"{organism}_{index_size}_create_index_star.stdout.log")
shell:
"(mkdir -p {params.output_dir}; \
chmod -R 777 {params.output_dir}; \
STAR \
--runMode genomeGenerate \
--sjdbOverhang {params.sjdbOverhang} \
--genomeDir {params.output_dir} \
--genomeFastaFiles {input.genome} \
--runThreadN {threads} \
--outFileNamePrefix {params.outFileNamePrefix} \
BIOPZ-Katsantoni Maria
committed
--sjdbGTFfile {input.gtf}) \
1> {log.stdout} 2> {log.stderr}"
BIOPZ-Gypas Foivos
committed
BIOPZ-Iborra de Toledo Paula
committed
rule extract_transcriptome:
"""
Create transcriptome from genome and gene annotations
"""
BIOPZ-Iborra de Toledo Paula
committed
input:
genome = lambda wildcards:
BIOPZ-Katsantoni Maria
committed
get_sample(
'genome',
search_id='organism',
search_value=wildcards.organism),
BIOPZ-Iborra de Toledo Paula
committed
gtf = lambda wildcards:
BIOPZ-Katsantoni Maria
committed
get_sample(
'gtf',
search_id='organism',
search_value=wildcards.organism)
BIOPZ-Iborra de Toledo Paula
committed
output:
transcriptome = os.path.join(
config['output_dir'],
"transcriptome",
"{organism}",
"transcriptome.fa")
BIOPZ-Iborra de Toledo Paula
committed
log:
stderr = os.path.join(
config['log_dir'],
"{organism}_extract_transcriptome.log"),
stdout = os.path.join(
config['log_dir'],
"{organism}_extract_transcriptome.log")
BIOPZ-Iborra de Toledo Paula
committed
singularity:
BIOPZ-Iborra de Toledo Paula
committed
shell:
"(gffread \
-w {output.transcriptome} \
-g {input.genome} {input.gtf}) \
1> {log.stdout} 2> {log.stderr}"
BIOPZ-Gypas Foivos
committed
rule concatenate_transcriptome_and_genome:
"""
Concatenate genome and transcriptome
"""
input:
transcriptome = os.path.join(
config['output_dir'],
"transcriptome",
"{organism}",
"transcriptome.fa"),
genome = lambda wildcards:
BIOPZ-Katsantoni Maria
committed
get_sample(
'genome',
search_id='organism',
search_value=wildcards.organism)
output:
genome_transcriptome = os.path.join(
config['output_dir'],
"transcriptome",
"{organism}",
"genome_transcriptome.fa")
singularity:
"docker://bash:5.0.16"
log:
stderr = os.path.join(
config['log_dir'],
"{organism}_concatenate_transcriptome_and_genome.stderr.log")
shell:
"(cat {input.transcriptome} {input.genome} \
1> {output.genome_transcriptome}) \
2> {log.stderr}"
rule create_index_salmon:
"""
Create index for Salmon quantification
"""
input:
config['output_dir'],
"transcriptome",
"{organism}",
"genome_transcriptome.fa"),
chr_names = lambda wildcards:
os.path.join(
config['star_indexes'],
BIOPZ-Katsantoni Maria
committed
get_sample('organism'),
get_sample("index_size"),
output:
index = directory(
os.path.join(
config['salmon_indexes'],
"{organism}",
"{kmer}",
BIOPZ-Katsantoni Maria
committed
"salmon.idx"))
BIOPZ-Katsantoni Maria
committed
kmerLen = "{kmer}"
"docker://zavolab/salmon:1.1.0-slim"
BIOPZ-Katsantoni Maria
committed
BIOPZ-Katsantoni Maria
committed
stderr = os.path.join(
BIOPZ-Katsantoni Maria
committed
"{organism}_{kmer}_create_index_salmon.stderr.log"),
stdout = os.path.join(
config['log_dir'],
"{organism}_{kmer}_create_index_salmon.stdout.log")
BIOPZ-Katsantoni Maria
committed
--transcripts {input.genome_transcriptome} \
--decoys {input.chr_names} \
--index {output.index} \
--kmerLen {params.kmerLen} \
BIOPZ-Katsantoni Maria
committed
--threads {threads}) \
1> {log.stdout} 2> {log.stderr}"
BIOPZ-Gypas Foivos
committed
rule create_index_kallisto:
BIOPZ-Katsantoni Maria
committed
"""
Create index for Kallisto quantification
"""
BIOPZ-Iborra de Toledo Paula
committed
transcriptome = os.path.join(
config['output_dir'],
"transcriptome",
"{organism}",
"transcriptome.fa")
output:
index = os.path.join(
config['kallisto_indexes'],
"{organism}",
BIOPZ-Katsantoni Maria
committed
"kallisto.idx")
params:
output_dir = os.path.join(
config['kallisto_indexes'],
BIOPZ-Katsantoni Maria
committed
"{organism}")
BIOPZ-Katsantoni Maria
committed
BIOPZ-Katsantoni Maria
committed
stderr = os.path.join(
BIOPZ-Katsantoni Maria
committed
"{organism}_create_index_kallisto.stderr.log"),
stdout = os.path.join(
config['log_dir'],
"{organism}_create_index_kallisto.stdout.log")
shell:
"(mkdir -p {params.output_dir}; \
chmod -R 777 {params.output_dir}; \
BIOPZ-Katsantoni Maria
committed
kallisto index -i {output.index} {input.transcriptome}) \
1> {log.stdout} 2> {log.stderr}"
BIOPZ-Gypas Foivos
committed
BIOPZ-Katsantoni Maria
committed
"""
Convert transcripts to BED12 format
"""
BIOPZ-Katsantoni Maria
committed
get_sample('gtf')
BIOPZ-Katsantoni Maria
committed
output:
bed12 = os.path.join(
config['output_dir'],
BIOPZ-Katsantoni Maria
committed
"full_transcripts_protein_coding.bed")
"docker://zavolab/zgtf:0.1"
BIOPZ-Katsantoni Maria
committed
BIOPZ-Katsantoni Maria
committed
stdout = os.path.join(
config['log_dir'],
"extract_transcripts_as_bed12.stdout.log"),
BIOPZ-Katsantoni Maria
committed
stderr = os.path.join(
BIOPZ-Katsantoni Maria
committed
"extract_transcripts_as_bed12.stderr.log")
"(gtf2bed12 \
--gtf {input.gtf} \
--transcript_type protein_coding \
--bed12 {output.bed12}); \
1> {log.stdout} 2> {log.stderr}"
rule index_genomic_alignment_samtools:
'''
Index genome bamfile using samtools
'''
input:
bam = os.path.join(
config["output_dir"],
"{sample}",
"map_genome",
"{sample}.{seqmode}.Aligned.sortedByCoord.out.bam"),
output:
bai = os.path.join(
config["output_dir"],
"{sample}",
"map_genome",
"{sample}.{seqmode}.Aligned.sortedByCoord.out.bam.bai")
singularity:
"docker://zavolab/samtools:1.10-slim"
threads: 1
log:
stderr = os.path.join(
config["log_dir"],
"{sample}",
"index_genomic_alignment_samtools.{seqmode}.stderr.log"),
stdout = os.path.join(
config["log_dir"],
"{sample}",
"index_genomic_alignment_samtools.{seqmode}.stdout.log")
shell:
"(samtools index {input.bam} {output.bai};) \
1> {log.stdout} 2> {log.stderr}"
BIOPZ-Katsantoni Maria
committed
"""
BIOPZ-Katsantoni Maria
committed
"""
bam = lambda wildcards:
expand(
os.path.join(
config['output_dir'],
"samples",
"{sample}",
"map_genome",
"{sample}.{seqmode}.Aligned.sortedByCoord.out.bam"),
sample=wildcards.sample,
BIOPZ-Katsantoni Maria
committed
seqmode=get_sample(
'seqmode',
search_id='index',
search_value=wildcards.sample)),
bai = lambda wildcards:
expand(
os.path.join(
config['output_dir'],
"samples",
"{sample}",
"map_genome",
"{sample}.{seqmode}.Aligned.sortedByCoord.out.bam.bai"),
sample=wildcards.sample,
BIOPZ-Katsantoni Maria
committed
seqmode=get_sample(
'seqmode',
search_id='index',
search_value=wildcards.sample)),
transcripts_bed12 = os.path.join(
config['output_dir'],
BIOPZ-Katsantoni Maria
committed
"full_transcripts_protein_coding.bed")
output:
TIN_score = os.path.join(
config['output_dir'],
BIOPZ-Katsantoni Maria
committed
"TIN_score.tsv")
BIOPZ-Katsantoni Maria
committed
sample = "{sample}"
BIOPZ-Katsantoni Maria
committed
stderr = os.path.join(
"samples",
"{sample}",
"calculate_TIN_scores.log")
BIOPZ-Katsantoni Maria
committed
BIOPZ-Katsantoni Maria
committed
"docker://zavolab/tin_score_calculation:0.2.0-slim"
BIOPZ-Katsantoni Maria
committed
BIOPZ-Katsantoni Maria
committed
"(tin_score_calculation.py \
-i {input.bam} \
-r {input.transcripts_bed12} \
-c 0 \
--names {params.sample} \
BIOPZ-Katsantoni Maria
committed
-n 100 > {output.TIN_score};) 2> {log.stderr}"
rule merge_TIN_scores:
"""
Merge TIN scores tables
"""
input:
TIN_score = expand(
os.path.join(
config['output_dir'],
"{sample}",
"TIN",
"TIN_score.tsv"),
BIOPZ-Katsantoni Maria
committed
sample=pd.unique(samples_table.index.values)),
output:
TIN_scores_merged = os.path.join(
config['output_dir'],
"TIN_scores_merged.tsv")
log:
stderr = os.path.join(
config['log_dir'],
"merge_TIN_scores.stderr.log"),
stdout = os.path.join(
config["log_dir"],
"merge_TIN_scores.stdout.log")
params:
TIN_score_merged_paths = ",".join(expand(
os.path.join(
config['output_dir'],
"{sample}",
"TIN",
"TIN_score.tsv"),
zip,
BIOPZ-Katsantoni Maria
committed
sample=[i for i in pd.unique(samples_table.index.values)],
seqmode=[get_sample('seqmode',
search_id='index',
search_value=i) for i in pd.unique(samples_table.index.values)]))
threads: 1
singularity:
"docker://zavolab/tin_score_calculation:0.2.0-slim"
shell:
"(tin_score_merge.py \
--input-files {params.TIN_score_merged_paths} \
--output-file {output.TIN_scores_merged}) \
1> {log.stdout} 2> {log.stderr}"
rule plot_TIN_scores:
"""
Generate TIN scores boxplots
"""
input:
TIN_scores_merged = os.path.join(
config['output_dir'],
"TIN_scores_merged.tsv"),
output:
TIN_boxplot_PNG = os.path.join(
config['output_dir'],
TIN_boxplot_PDF = os.path.join(
config['output_dir'],
params:
TIN_boxplot_prefix = os.path.join(
config['output_dir'],
log:
stderr = os.path.join(
config['log_dir'],
"plot_TIN_scores.stderr.log"),
stdout = os.path.join(
config["log_dir"],
"plot_TIN_scores.stdout.log")
threads: 1
singularity:
"docker://zavolab/tin_score_calculation:0.2.0-slim"
shell:
"(tin_score_plot.py \
--input-file {input.TIN_scores_merged} \
--output-file-prefix {params.TIN_boxplot_prefix}) \
1> {log.stdout} 2> {log.stderr}"
rule salmon_quantmerge_genes:
BIOPZ-Katsantoni Maria
committed
'''
Merge gene quantifications into a single file
'''
BIOPZ-Katsantoni Maria
committed
salmon_in = expand(
os.path.join(
config["output_dir"],
BIOPZ-Katsantoni Maria
committed
"{sample}",
"{sample}.salmon.{seqmode}",
"quant.sf"),
BIOPZ-Katsantoni Maria
committed
sample=pd.unique(samples_table.index.values),
seqmode=[get_sample(
'seqmode',
search_id='index',
search_value=i)
for i in pd.unique(samples_table.index.values)])
BIOPZ-Katsantoni Maria
committed
BIOPZ-Katsantoni Maria
committed
salmon_out = os.path.join(
BIOPZ-Katsantoni Maria
committed
"summary_salmon",
"quantmerge",
"genes_{salmon_merge_on}.tsv")
params:
BIOPZ-Katsantoni Maria
committed
os.path.join(
config["output_dir"],
BIOPZ-Katsantoni Maria
committed
"{sample}",
BIOPZ-Katsantoni Maria
committed
sample=[i for i in pd.unique(samples_table.index.values)],
seqmode=[get_sample(
'seqmode',
search_id='index',
search_value=i)
for i in pd.unique(samples_table.index.values)]),
BIOPZ-Katsantoni Maria
committed
sample_name_list = expand(
"{sample}",
BIOPZ-Katsantoni Maria
committed
sample=pd.unique(samples_table.index.values)),
salmon_merge_on = "{salmon_merge_on}"
BIOPZ-Katsantoni Maria
committed
BIOPZ-Katsantoni Maria
committed
stderr = os.path.join(
config["log_dir"],
"salmon_quantmerge_genes_{salmon_merge_on}.stderr.log"),
stdout = os.path.join(
config["log_dir"],
"salmon_quantmerge_genes_{salmon_merge_on}.stdout.log")
threads: 1
singularity:
"docker://zavolab/salmon:1.1.0-slim"
BIOPZ-Katsantoni Maria
committed
shell:
"(salmon quantmerge \
--genes \
--names {params.sample_name_list} \
--column {params.salmon_merge_on} \
BIOPZ-Katsantoni Maria
committed
--output {output.salmon_out};) \
1> {log.stdout} 2> {log.stderr}"
rule salmon_quantmerge_transcripts:
BIOPZ-Katsantoni Maria
committed
'''
Merge gene quantifications into a single file
'''
BIOPZ-Katsantoni Maria
committed
salmon_in = expand(
os.path.join(
config["output_dir"],
BIOPZ-Katsantoni Maria
committed
"{sample}",
BIOPZ-Katsantoni Maria
committed
"quant.sf"),
BIOPZ-Katsantoni Maria
committed
sample=[i for i in pd.unique(samples_table.index.values)],
seqmode=[get_sample(
'seqmode',
search_id='index',
search_value=i)
for i in pd.unique(samples_table.index.values)])
BIOPZ-Katsantoni Maria
committed
BIOPZ-Katsantoni Maria
committed
salmon_out = os.path.join(
BIOPZ-Katsantoni Maria
committed
"summary_salmon",
"quantmerge",
"transcripts_{salmon_merge_on}.tsv")
params:
BIOPZ-Katsantoni Maria
committed
os.path.join(
config["output_dir"],
BIOPZ-Katsantoni Maria
committed
"{sample}",
BIOPZ-Katsantoni Maria
committed
sample=[i for i in pd.unique(samples_table.index.values)],
seqmode=[get_sample(
'seqmode',
search_id='index',
search_value=i)
for i in pd.unique(samples_table.index.values)]),
BIOPZ-Katsantoni Maria
committed
sample_name_list = expand(
"{sample}",
BIOPZ-Katsantoni Maria
committed
sample=pd.unique(samples_table.index.values)),
salmon_merge_on = "{salmon_merge_on}"
BIOPZ-Katsantoni Maria
committed
BIOPZ-Katsantoni Maria
committed
stderr = os.path.join(
config["log_dir"],
"salmon_quantmerge_transcripts_{salmon_merge_on}.stderr.log"),
stdout = os.path.join(
config["log_dir"],
"salmon_quantmerge_transcripts_{salmon_merge_on}.stdout.log")
threads: 1
singularity:
"docker://zavolab/salmon:1.1.0-slim"
BIOPZ-Katsantoni Maria
committed
shell:
"(salmon quantmerge \
--names {params.sample_name_list} \
--column {params.salmon_merge_on} \
BIOPZ-Katsantoni Maria
committed
--output {output.salmon_out}) \
rule star_rpm:
'''
Create stranded bedgraph coverage with STARs RPM normalisation
'''
input:
bam = lambda wildcards:
expand(
os.path.join(
config["output_dir"],
"samples",
"{sample}",
"map_genome",
"{sample}.{seqmode}.Aligned.sortedByCoord.out.bam"),
sample=wildcards.sample,
BIOPZ-Katsantoni Maria
committed
seqmode=get_sample(
'seqmode',
search_id='index',
search_value=wildcards.sample)),
bai = lambda wildcards:
expand(
os.path.join(
config["output_dir"],
"samples",
"{sample}",
"map_genome",
"{sample}.{seqmode}.Aligned.sortedByCoord.out.bam.bai"),
sample=wildcards.sample,
BIOPZ-Katsantoni Maria
committed
seqmode=get_sample(
'seqmode',
search_id='index',
search_value=wildcards.sample))
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
output:
str1 = os.path.join(
config["output_dir"],
"samples",
"{sample}",
"STAR_coverage",
"{sample}_Signal.Unique.str1.out.bg"),
str2 = os.path.join(
config["output_dir"],
"samples",
"{sample}",
"STAR_coverage",
"{sample}_Signal.UniqueMultiple.str1.out.bg"),
str3 = os.path.join(
config["output_dir"],
"samples",
"{sample}",
"STAR_coverage",
"{sample}_Signal.Unique.str2.out.bg"),
str4 = os.path.join(
config["output_dir"],
"samples",
"{sample}",
"STAR_coverage",
"{sample}_Signal.UniqueMultiple.str2.out.bg")
params:
out_dir = lambda wildcards, output:
os.path.dirname(output.str1),
prefix = lambda wildcards, output:
os.path.join(
os.path.dirname(output.str1),
str(wildcards.sample) + "_"),
stranded = "Stranded"
singularity:
"docker://zavolab/star:2.7.3a-slim"
log:
stderr = os.path.join(
config["log_dir"],
"samples",
"{sample}",
"star_rpm.stderr.log"),
stdout = os.path.join(
config["log_dir"],
"samples",
"{sample}",
"star_rpm.stdout.log")
threads: 4
shell:
"(mkdir -p {params.out_dir}; \
chmod -R 777 {params.out_dir}; \
STAR \
--runMode inputAlignmentsFromBAM \
--runThreadN {threads} \
--inputBAMfile {input.bam} \
--outWigType bedGraph \
--outWigStrand {params.stranded} \
--outWigNorm RPM \
--outFileNamePrefix {params.prefix}) \
1> {log.stdout} 2> {log.stderr}"
rule rename_star_rpm_for_alfa:
input:
plus = lambda wildcards:
expand(
os.path.join(
config["output_dir"],
"samples",
"{sample}",
"STAR_coverage",
"{sample}_Signal.{unique}.{plus}.out.bg"),
sample=wildcards.sample,
unique=wildcards.unique,
BIOPZ-Katsantoni Maria
committed
plus=get_sample(
'alfa_plus',
search_id='index',
search_value=wildcards.sample)),
minus = lambda wildcards:
expand(
os.path.join(
config["output_dir"],
"samples",
"{sample}",
"STAR_coverage",
"{sample}_Signal.{unique}.{minus}.out.bg"),
sample=wildcards.sample,
unique=wildcards.unique,
BIOPZ-Katsantoni Maria
committed
minus=get_sample(
'alfa_minus',
search_id='index',
search_value=wildcards.sample))
output:
plus = os.path.join(
config["output_dir"],
"samples",
"{sample}",
"ALFA",
"{unique}",
"{sample}.{unique}.plus.bg"),
minus = os.path.join(
config["output_dir"],
"samples",
"{sample}",
"ALFA",
"{unique}",
"{sample}.{unique}.minus.bg")
params:
orientation = lambda wildcards:
BIOPZ-Katsantoni Maria
committed
get_sample(
'kallisto_directionality',
search_id='index',
search_value=wildcards.sample),
log:
stderr = os.path.join(
config["log_dir"],
"samples",
"{sample}",
"rename_star_rpm_for_alfa__{unique}.stderr.log"),
stdout = os.path.join(
config["log_dir"],
"samples",
"{sample}",
"rename_star_rpm_for_alfa__{unique}.stdout.log")
singularity:
"docker://bash:5.0.16"
shell:
"(cp {input.plus} {output.plus}; \
cp {input.minus} {output.minus};) \
1>{log.stdout} 2>{log.stderr}"
rule generate_alfa_index:
''' Generate ALFA index files from sorted GTF file '''
input:
BIOPZ-Katsantoni Maria
committed
get_sample(
'gtf',
search_id='organism',
search_value=wildcards.organism),
chr_len = os.path.join(
config["star_indexes"],
"{organism}",
"{index_size}",
"STAR_index",
"chrNameLength.txt"),
output:
index_stranded = os.path.join(
config["alfa_indexes"],
"{organism}",
"{index_size}",
"ALFA",
index_unstranded = os.path.join(
config["alfa_indexes"],
"{organism}",
"{index_size}",
"ALFA",
"sorted_genes.unstranded.ALFA_index")
params:
genome_index = "sorted_genes",
out_dir = lambda wildcards, output:
os.path.dirname(output.index_stranded)
log:
os.path.join(
config["log_dir"],
"{organism}_{index_size}_generate_alfa_index.log")
"(alfa -a {input.gtf} \
-g {params.genome_index} \
--chr_len {input.chr_len} \
-p {threads} \
-o {params.out_dir}) &> {log}"
'''
Run ALFA from stranded bedgraph files
'''
input:
plus = os.path.join(
config["output_dir"],
gtf = lambda wildcards:
os.path.join(
config["alfa_indexes"],
BIOPZ-Katsantoni Maria
committed
get_sample(
'organism',
search_id='index',
search_value=wildcards.sample),
get_sample(
'index_size',
search_id='index',
search_value=wildcards.sample),
"ALFA",
"sorted_genes.stranded.ALFA_index")
output:
biotypes = os.path.join(
config["output_dir"],
"ALFA_plots.Biotypes.pdf"),
categories = os.path.join(
config["output_dir"],
"ALFA_plots.Categories.pdf"),
table = os.path.join(
config["output_dir"],
out_dir = lambda wildcards, output:
os.path.dirname(output.biotypes),
BIOPZ-Katsantoni Maria
committed
get_sample(
'alfa_directionality',
search_id='index',
search_value=wildcards.sample),
genome_index = lambda wildcards, input:
os.path.abspath(
os.path.join(
os.path.dirname(input.gtf),
"sorted_genes")),
"samples",
"{sample}",
"alfa_qc.{unique}.log")
"(cd {params.out_dir}; \
alfa \
-g {params.genome_index} \
--bedgraph {params.plus} {params.minus} {params.name} \
'''
Run ALFA from stranded bedgraph files on all samples
'''
tables = lambda wildcards:
expand(
os.path.join(
config["output_dir"],
"samples",
"{sample}",
"ALFA",
"{unique}",
BIOPZ-Katsantoni Maria
committed
sample=pd.unique(samples_table.index.values),
output:
biotypes = os.path.join(
config["output_dir"],
"ALFA",
"ALFA_plots.Biotypes.pdf"),
categories = os.path.join(
config["output_dir"],
"ALFA",
out_dir = lambda wildcards, output:
os.path.dirname(output.biotypes)
log:
os.path.join(
config["log_dir"],
"alfa_qc_all_samples.{unique}.log")
"(alfa -c {input.tables} -o {params.out_dir}) &> {log}"
expand(
os.path.join(
config["output_dir"],
"ALFA",
"{unique}",
"ALFA_plots.{annotation}.pdf"),
unique=["Unique", "UniqueMultiple"],
annotation=["Categories", "Biotypes"])
log:
os.path.join(
config["log_dir"],
"alfa_qc_all_samples.concat.log")
"docker://zavolab/imagemagick:7.0.8"
"(convert -append -density {params.density} \
{input} {output}) &> {log}"
rule prepare_multiqc_config:
'''
Prepare config for the MultiQC
'''
input:
script = os.path.join(
workflow.basedir,
"workflow",
"scripts",
output:
multiqc_config = os.path.join(
config["output_dir"],
"multiqc_config.yaml")
logo_path = config['report_logo'],
multiqc_intro_text = config['report_description'],
url = config['report_url']
log:
stderr = os.path.join(
config["log_dir"],
"prepare_multiqc_config.stderr.log"),
stdout = os.path.join(
config["log_dir"],
"prepare_multiqc_config.stdout.log")
shell:
--intro-text '{params.multiqc_intro_text}' \
--custom-logo {params.logo_path} \
--url '{params.url}') \
BIOPZ-Katsantoni Maria
committed
sample=pd.unique(samples_table.index.values),
mate="fq1"),
fastqc_pe = expand(
os.path.join(
config['output_dir'],
"samples",
"{sample}",
"fastqc",
"{mate}"),
BIOPZ-Katsantoni Maria
committed
sample=[i for i in pd.unique(
samples_table[samples_table['seqmode'] == 'pe'].index.values)],
mate="fq2"),
pseudoalignment = expand(
os.path.join(
config['output_dir'],
BIOPZ-Katsantoni Maria
committed
sample=[i for i in pd.unique(samples_table.index.values)],
seqmode=[get_sample('seqmode', search_id='index', search_value=i)
for i in pd.unique(samples_table.index.values)]),
config["output_dir"],
"multiqc_config.yaml")
output:
multiqc_report = directory(
os.path.join(
config["output_dir"],
"multiqc_summary"))
results_dir = os.path.join(
config["output_dir"]),
log_dir = config["log_dir"]
stderr = os.path.join(
config["log_dir"],
stdout = os.path.join(
config["log_dir"],
singularity:
"docker://ewels/multiqc:1.7"
shell:
"(multiqc \
--outdir {output.multiqc_report} \
--config {input.multiqc_config} \
{params.results_dir} \
{params.log_dir};) \
1> {log.stdout} 2> {log.stderr}"
sort bedGraphs in order to work with bedGraphtobigWig
config["output_dir"],
"samples",
"{sample}",
"ALFA",
"{unique}",
sorted_bg = os.path.join(
config["output_dir"],
"samples",
"{sample}",
"bigWig",
"{unique}",
"{sample}_{unique}_{strand}.sorted.bg")
singularity:
"docker://cjh4zavolab/bedtools:2.27"
stderr = os.path.join(
config["log_dir"],
"sort_bg_{unique}_{strand}.stderr.log")
"(sortBed \
-i {input.bg} \
> {output.sorted_bg};) 2> {log.stderr}"
bedGraphtobigWig, for viewing in genome browsers
sorted_bg = os.path.join(
config["output_dir"],
"samples",
"{sample}",
"bigWig",
"{unique}",
"{sample}_{unique}_{strand}.sorted.bg"),
chr_sizes = lambda wildcards:
os.path.join(
config['star_indexes'],
BIOPZ-Katsantoni Maria
committed
get_sample(
'organism',
search_id='index',
search_value=wildcards.sample),
get_sample(
'index_size',
search_id='index',
search_value=wildcards.sample),
"STAR_index",
"chrNameLength.txt")
bigWig = os.path.join(
config["output_dir"],
"samples",
"{sample}",
"bigWig",
"{unique}",
"{sample}_{unique}_{strand}.bw")
singularity:
"docker://zavolab/bedgraphtobigwig:4-slim"
stderr = os.path.join(
config["log_dir"],
"samples",
"{sample}",
"bigwig_{unique}_{strand}.stderr.log"),
stdout = os.path.join(
config["log_dir"],
"samples",
"{sample}",
"bigwig_{unique}_{strand}.stdout.log")
"(bedGraphToBigWig \
{input.sorted_bg} \
{input.chr_sizes} \
{output.bigWig};) \
1> {log.stdout} 2> {log.stderr}"