Merge branch 'alfa' into 'master'

Alfa integration

See merge request zavolan_group/pipelines/rnaseqpipeline!41

.gitignore

@@ -332,4 +332,5 @@ runs/.*
 .snakemake/
 logs/
 results/
+!tests/test_alfa/results/
 .java/

.gitlab-ci.yml

@@ -14,6 +14,7 @@ test:
     # add script tests here
     - bash tests/test_scripts_labkey_to_snakemake_table/test.sh
     - bash tests/test_scripts_labkey_to_snakemake_api/test.sh
+    - bash tests/test_alfa/test.sh
     # add integration tests here
     - bash tests/test_create_dag_image/test.sh
     - bash tests/test_create_rule_graph/test.sh

Snakefile

@@ -4,6 +4,7 @@ import os
 import sys
 
 import pandas as pd
+import shutil
 
 # Get sample table
 samples_table = pd.read_csv(
@@ -16,7 +17,7 @@ samples_table = pd.read_csv(
 )
 
 # Global config
-localrules: finish
+localrules: finish, rename_star_rpm_for_alfa
 
 # Create log directories
 os.makedirs(
@@ -99,7 +100,23 @@ rule finish:
                 zip,
                 sample=[i for i in list(samples_table.index.values)],
                 seqmode=[samples_table.loc[i, 'seqmode']
-                        for i in list(samples_table.index.values)])
+                        for i in list(samples_table.index.values)]),
+        alfa_reports = expand(os.path.join(
+            config["output_dir"],
+            "{seqmode}",
+            "{sample}",
+            "ALFA",
+            "ALFA_plots.Biotypes.pdf"),
+            zip,
+            sample= [i for i in list(samples_table.index.values)],
+            seqmode= [
+                samples_table.loc[i,"seqmode"] 
+                for i in list(samples_table.index.values)]),
+        alfa_all_samples = os.path.join(
+            config["output_dir"],
+            "ALFA",
+            "ALFA_plots.Categories.pdf")
+
 
 
 
@@ -441,6 +458,47 @@ rule star_rpm:
         """
 
 
+rule rename_star_rpm_for_alfa:
+    input:
+        str1 = os.path.join(
+            config["output_dir"],
+            "{seqmode}",
+            "{sample}",
+            "STAR_coverage",
+            "{sample}_Signal.UniqueMultiple.str1.out.bg"),
+        str2 = os.path.join(
+            config["output_dir"],
+            "{seqmode}",
+            "{sample}",
+            "STAR_coverage",
+            "{sample}_Signal.UniqueMultiple.str2.out.bg")
+    
+    output:
+        plus = os.path.join(
+            config["output_dir"],
+            "{seqmode}",
+            "{sample}",
+            "ALFA",
+            "{sample}_Signal.UniqueMultiple.out.plus.bg"),
+        minus = os.path.join(
+            config["output_dir"],
+            "{seqmode}",
+            "{sample}",
+            "ALFA",
+            "{sample}_Signal.UniqueMultiple.out.minus.bg")
+    
+    params:
+        orientation = lambda wildcards: samples_table.loc[wildcards.sample, "kallisto_directionality"]
+    
+    run:
+        if params['orientation'] == "--fr":
+            shutil.copy2(input['str1'], output['plus'])
+            shutil.copy2(input['str2'], output['minus'])
+        elif params['orientation'] == "--rf":
+            shutil.copy2(input['str1'], output['minus'])
+            shutil.copy2(input['str2'], output['plus'])
+
+
 rule calculate_TIN_scores:
     """
         Caluclate transcript integrity (TIN) score
@@ -612,4 +670,162 @@ rule salmon_quantmerge_transcripts:
         --names {params.sample_name_list} \
         --column {params.salmon_merge_on} \
         --output {output.salmon_out}) \
-        1> {log.stdout} 2> {log.stderr}"
+        1> {log.stdout} 2> {log.stderr}"
\ No newline at end of file
+
+
+#################################################################################
+### ALFA: Annotation Landscape For Aligned reads
+#################################################################################
+
+directionality = {"--fr": "fr-firststrand", "--rf": "fr-secondstrand"}
+
+
+rule generate_alfa_index:
+    ''' Generate ALFA index files from sorted GTF file '''
+    input:
+        gtf = lambda wildcards: samples_table["gtf"][samples_table["organism"]==wildcards.organism][0],
+        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", 
+            "sorted_genes.stranded.ALFA_index"),
+        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)
+
+    threads: 4
+
+    singularity: 
+        "docker://zavolab/alfa:1.1.1"
+
+    log: 
+        os.path.join(config["log_dir"], "{organism}_{index_size}_generate_alfa_index.log")
+
+    shell:
+        """
+        alfa -a {input.gtf} \
+            -g {params.genome_index} \
+            --chr_len {input.chr_len} \
+            -p {threads} \
+            -o {params.out_dir} &> {log}
+        """
+
+
+rule alfa_qc:
+    ''' Run ALFA from stranded bedgraph files '''
+    input:
+        plus = os.path.join(
+            config["output_dir"],
+            "{seqmode}",
+            "{sample}",
+            "ALFA",
+            "{sample}_Signal.UniqueMultiple.out.plus.bg"),
+        minus = os.path.join(
+            config["output_dir"],
+            "{seqmode}",
+            "{sample}",
+            "ALFA",
+            "{sample}_Signal.UniqueMultiple.out.minus.bg"),
+        gtf = lambda wildcards: os.path.join(config["alfa_indexes"], 
+            samples_table.loc[wildcards.sample, "organism"], 
+            str(samples_table.loc[wildcards.sample, "index_size"]), 
+            "ALFA", 
+            "sorted_genes.stranded.ALFA_index")
+
+    output:
+        biotypes = os.path.join(
+            config["output_dir"],
+            "{seqmode}",
+            "{sample}",
+            "ALFA",
+            "ALFA_plots.Biotypes.pdf"),
+        categories = os.path.join(
+            config["output_dir"],
+            "{seqmode}",
+            "{sample}",
+            "ALFA",
+            "ALFA_plots.Categories.pdf"),
+        table = os.path.join(
+            config["output_dir"],
+            "{seqmode}",
+            "{sample}",
+            "ALFA",
+            "{sample}.ALFA_feature_counts.tsv")
+
+    params:
+        out_dir = lambda wildcards, output: os.path.dirname(output.biotypes),
+        alfa_orientation = lambda wildcards: directionality[samples_table.loc[wildcards.sample, "kallisto_directionality"]],
+        in_file_plus = lambda wildcards, input: os.path.basename(input.plus),
+        in_file_minus = lambda wildcards, input: os.path.basename(input.minus),
+        genome_index = lambda wildcards, input: os.path.abspath(os.path.join(os.path.dirname(input.gtf), "sorted_genes")),
+        name = "{sample}"
+
+    singularity:
+        "docker://zavolab/alfa:1.1.1"
+
+    log: 
+        os.path.abspath(os.path.join(
+            config["log_dir"], 
+            "{seqmode}", 
+            "{sample}", 
+            "alfa_qc.log"))
+
+    shell:
+        """ 
+        cd {params.out_dir}; \
+        (alfa -g {params.genome_index} \
+            --bedgraph {params.in_file_plus} {params.in_file_minus} {params.name} \
+            -s {params.alfa_orientation}) &> {log}
+        """
+
+
+rule alfa_qc_all_samples:
+    ''' Run ALFA from stranded bedgraph files on all samples '''
+    input:
+        tables = [os.path.join(
+            config["output_dir"],
+            samples_table.loc[sample1, "seqmode"],
+            str(sample1),
+            "ALFA",
+            sample1 + ".ALFA_feature_counts.tsv")
+            for sample1 in list(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",
+            "ALFA_plots.Categories.pdf")
+
+    params:
+        out_dir = lambda wildcards, output: os.path.dirname(output.biotypes)
+
+    log: 
+        os.path.abspath(
+            os.path.join(config["log_dir"], 
+            "alfa_qc_all_samples.log"))
+
+    singularity:
+        "docker://zavolab/alfa:1.1.1"
+
+    shell:
+        """
+        (alfa -c {input.tables} -o {params.out_dir}) &> {log}
+        """

scripts/labkey_to_snakemake.py

@@ -275,6 +275,7 @@ def main():
   kallisto_indexes: "results/kallisto_indexes/"
   salmon_indexes: "results/salmon_indexes/"
   star_indexes: "results/star_indexes/"
+  alfa_indexes: "results/alfa_indexes/"
 ...''')
 
     sys.stdout.write('Create snakemake table finished successfully...\n')

tests/input_files/config.yaml

@@ -5,4 +5,5 @@
   kallisto_indexes: "results/kallisto_indexes/"
   salmon_indexes: "results/salmon_indexes/"
   star_indexes: "results/star_indexes/"
+  alfa_indexes: "results/alfa_indexes/"
 ...

tests/input_files/config_alfa.yaml

+---
+  samples: "../input_files/samples_alfa.tsv"
+  output_dir: "results/"
+  log_dir: "logs/"
+  kallisto_indexes: "results/kallisto_indexes/"
+  salmon_indexes: "results/salmon_indexes/"
+  star_indexes: "results/star_indexes/"
+  alfa_indexes: "results/alfa_indexes/"
+...

tests/input_files/homo_sapiens/annotation.gtf

@@ -15,8 +15,8 @@
 1-10000-20000	havana	exon	2976	3053	.	+	.	gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "4"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001799933"; exon_version "2"; tag "basic"; transcript_support_level "NA";
 1-10000-20000	havana	exon	3222	3375	.	+	.	gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "5"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001746346"; exon_version "2"; tag "basic"; transcript_support_level "NA";
 1-10000-20000	havana	exon	3454	3671	.	+	.	gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "6"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001863096"; exon_version "1"; tag "basic"; transcript_support_level "NA";
-1-10000-20000	havana	gene	4405	8367	.	-	.	gene_id	"ENSG00000227232"; gene_version	"5"; gene_name "WASH7P"; gene_source "havana"; gene_biotype "unprocessed_pseudogene";
-1-10000-20000	havana	transcript	4405	8367	.	-	.	gene_id "ENSG00000227232"; gene_version "5"; transcript_id "ENST00000488147"; transcript_version "1"; gene_name	"WASH7P"; gene_source "havana"; gene_biotype "unprocessed_pseudogene"; transcript_name "WASH7P-201"; transcript_source "havana"; transcript_biotype "unprocessed_pseudogene"; tag "basic"; transcript_support_level "NA";
+1-10000-20000	havana	gene	4405	8367	.	-	.	gene_id "ENSG00000227232"; gene_version "5"; gene_name "WASH7P"; gene_source "havana"; gene_biotype "unprocessed_pseudogene";
+1-10000-20000	havana	transcript	4405	8367	.	-	.	gene_id "ENSG00000227232"; gene_version "5"; transcript_id "ENST00000488147"; transcript_version "1"; gene_name "WASH7P"; gene_source "havana"; gene_biotype "unprocessed_pseudogene"; transcript_name "WASH7P-201"; transcript_source "havana"; transcript_biotype "unprocessed_pseudogene"; tag "basic"; transcript_support_level "NA";
 1-10000-20000	havana	exon	8269	8367	.	-	.	gene_id "ENSG00000227232"; gene_version "5"; transcript_id "ENST00000488147"; transcript_version "1"; exon_number "3"; gene_name "WASH7P"; gene_source "havana"; gene_biotype "unprocessed_pseudogene"; transcript_name "WASH7P-201"; transcript_source "havana"; transcript_biotype "unprocessed_pseudogene"; exon_id "ENSE00003477500"; exon_version "1"; tag "basic"; transcript_support_level "NA";
 1-10000-20000	havana	exon	7916	8062	.	-	.	gene_id "ENSG00000227232"; gene_version "5"; transcript_id "ENST00000488147"; transcript_version "1"; exon_number "4"; gene_name "WASH7P"; gene_source "havana"; gene_biotype "unprocessed_pseudogene"; transcript_name "WASH7P-201"; transcript_source "havana"; transcript_biotype "unprocessed_pseudogene"; exon_id "ENSE00003565697"; exon_version "1"; tag "basic"; transcript_support_level "NA";
 1-10000-20000	havana	exon	7607	7743	.	-	.	gene_id "ENSG00000227232"; gene_version "5"; transcript_id "ENST00000488147"; transcript_version "1"; exon_number "5"; gene_name "WASH7P"; gene_source "havana"; gene_biotype "unprocessed_pseudogene"; transcript_name "WASH7P-201"; transcript_source "havana"; transcript_biotype "unprocessed_pseudogene"; exon_id "ENSE00003475637"; exon_version "1"; tag "basic"; transcript_support_level "NA";

tests/input_files/homo_sapiens/quick_start.gtf

+Chr1	ensembl_havana	gene	251	1250	.	+	.	gene_id "ENSMUSG00000051951"; gene_biotype "protein_coding";
+Chr1	ensembl_havana	transcript	251	1250	.	+	.	gene_id "ENSMUSG00000051951"; gene_biotype "protein_coding";
+Chr1	ensembl_havana	exon	376	1000	.	+	.	gene_id "ENSMUSG00000051951"; gene_biotype "protein_coding";
+Chr1	ensembl_havana	CDS	376	1000	.	+	0	gene_id "ENSMUSG00000051951"; gene_biotype "protein_coding";
+Chr1	ensembl_havana	five_prime_utr	251	375	.	+	.	gene_id "ENSMUSG00000051951"; gene_biotype "protein_coding";
+Chr1	ensembl_havana	three_prime_utr	1001	1250	.	+	.	gene_id "ENSMUSG00000051951"; gene_biotype "protein_coding";

tests/input_files/samples_alfa.tsv

+sample	seqmode	fq1	index_size	kmer	fq2	fq1_3p	fq1_5p	fq2_3p	fq2_5p	organism	gtf	gtf_filtered	genome	tr_fasta_filtered	sd	mean	multimappers	soft_clip	pass_mode	libtype	kallisto_directionality	fq1_polya	fq2_polya
+paired_end_R1_on_plus_sense	paired_end	XXXXXXXXXXXXX	75	31	XXXXXXXXXXXXX	GATCGGAAGAGCACA	XXXXXXXXXXXXX	AGATCGGAAGAGCGT	XXXXXXXXXXXXX	homo_sapiens	../input_files/homo_sapiens/quick_start.gtf	../input_files/homo_sapiens/quick_start.gtf	XXXXXXXXXXXXX	XXXXXXXXXXXXX	100	250	10	EndToEnd	None	A	--fr	AAAAAAAAAAAAAAAAA	TTTTTTTTTTTTTTTTT
+paired_end_R1_on_plus_antisense	paired_end	XXXXXXXXXXXXX	75	31	XXXXXXXXXXXXX	AGATCGGAAGAGCACA	XXXXXXXXXXXXX	AGATCGGAAGAGCGT	XXXXXXXXXXXXX	homo_sapiens	../input_files/homo_sapiens/quick_start.gtf	../input_files/homo_sapiens/quick_start.gtf	XXXXXXXXXXXXX	XXXXXXXXXXXXX	100	250	10	EndToEnd	None	A	--rf	AAAAAAAAAAAAAAAAA	TTTTTTTTTTTTTTTTT
+paired_end_R1_on_minus_sense	paired_end	XXXXXXXXXXXXX	75	31	XXXXXXXXXXXXX	AGATCGGAAGAGCACA	XXXXXXXXXXXXX	AGATCGGAAGAGCGT	XXXXXXXXXXXXX	homo_sapiens	../input_files/homo_sapiens/quick_start.gtf	../input_files/homo_sapiens/quick_start.gtf	XXXXXXXXXXXXX	XXXXXXXXXXXXX	100	250	10	EndToEnd	None	A	--fr	AAAAAAAAAAAAAAAAA	TTTTTTTTTTTTTTTTT
+paired_end_R1_on_minus_antisense	paired_end	XXXXXXXXXXXXX	75	31	XXXXXXXXXXXXX	AGATCGGAAGAGCACA	XXXXXXXXXXXXX	AGATCGGAAGAGCGT	XXXXXXXXXXXXX	homo_sapiens	../input_files/homo_sapiens/quick_start.gtf	../input_files/homo_sapiens/quick_start.gtf	XXXXXXXXXXXXX	XXXXXXXXXXXXX	100	250	10	EndToEnd	None	A	--rf	AAAAAAAAAAAAAAAAA	TTTTTTTTTTTTTTTTT

tests/test_alfa/expected_output.md5

+90e42aa46890e9cd0a47800428699fbf  results/alfa_indexes/homo_sapiens/75/ALFA/sorted_genes.stranded.ALFA_index
+ccc3cf5a57fddb0d469e597d4376b1bf  results/alfa_indexes/homo_sapiens/75/ALFA/sorted_genes.unstranded.ALFA_index
+5e90c760710980f4f4866dbe9aa32c6c  results/paired_end/paired_end_R1_on_minus_antisense/STAR_coverage/paired_end_R1_on_minus_antisense_Signal.UniqueMultiple.str1.out.bg
+8e23d52d7f635d927e292174f33168eb  results/paired_end/paired_end_R1_on_minus_antisense/STAR_coverage/paired_end_R1_on_minus_antisense_Signal.UniqueMultiple.str2.out.bg
+5e90c760710980f4f4866dbe9aa32c6c  results/paired_end/paired_end_R1_on_minus_antisense/STAR_coverage/paired_end_R1_on_minus_antisense_Signal.Unique.str1.out.bg
+8e23d52d7f635d927e292174f33168eb  results/paired_end/paired_end_R1_on_minus_antisense/STAR_coverage/paired_end_R1_on_minus_antisense_Signal.Unique.str2.out.bg
+5e90c760710980f4f4866dbe9aa32c6c  results/paired_end/paired_end_R1_on_minus_sense/STAR_coverage/paired_end_R1_on_minus_sense_Signal.UniqueMultiple.str1.out.bg
+8e23d52d7f635d927e292174f33168eb  results/paired_end/paired_end_R1_on_minus_sense/STAR_coverage/paired_end_R1_on_minus_sense_Signal.UniqueMultiple.str2.out.bg
+5e90c760710980f4f4866dbe9aa32c6c  results/paired_end/paired_end_R1_on_minus_sense/STAR_coverage/paired_end_R1_on_minus_sense_Signal.Unique.str1.out.bg
+8e23d52d7f635d927e292174f33168eb  results/paired_end/paired_end_R1_on_minus_sense/STAR_coverage/paired_end_R1_on_minus_sense_Signal.Unique.str2.out.bg
+8e23d52d7f635d927e292174f33168eb  results/paired_end/paired_end_R1_on_plus_antisense/STAR_coverage/paired_end_R1_on_plus_antisense_Signal.UniqueMultiple.str1.out.bg
+5e90c760710980f4f4866dbe9aa32c6c  results/paired_end/paired_end_R1_on_plus_antisense/STAR_coverage/paired_end_R1_on_plus_antisense_Signal.UniqueMultiple.str2.out.bg
+8e23d52d7f635d927e292174f33168eb  results/paired_end/paired_end_R1_on_plus_antisense/STAR_coverage/paired_end_R1_on_plus_antisense_Signal.Unique.str1.out.bg
+5e90c760710980f4f4866dbe9aa32c6c  results/paired_end/paired_end_R1_on_plus_antisense/STAR_coverage/paired_end_R1_on_plus_antisense_Signal.Unique.str2.out.bg
+8e23d52d7f635d927e292174f33168eb  results/paired_end/paired_end_R1_on_plus_sense/STAR_coverage/paired_end_R1_on_plus_sense_Signal.UniqueMultiple.str1.out.bg
+5e90c760710980f4f4866dbe9aa32c6c  results/paired_end/paired_end_R1_on_plus_sense/STAR_coverage/paired_end_R1_on_plus_sense_Signal.UniqueMultiple.str2.out.bg
+8e23d52d7f635d927e292174f33168eb  results/paired_end/paired_end_R1_on_plus_sense/STAR_coverage/paired_end_R1_on_plus_sense_Signal.Unique.str1.out.bg
+5e90c760710980f4f4866dbe9aa32c6c  results/paired_end/paired_end_R1_on_plus_sense/STAR_coverage/paired_end_R1_on_plus_sense_Signal.Unique.str2.out.bg
+5e90c760710980f4f4866dbe9aa32c6c  results/paired_end/paired_end_R1_on_minus_antisense/ALFA/paired_end_R1_on_minus_antisense_Signal.UniqueMultiple.out.minus.bg
+8e23d52d7f635d927e292174f33168eb  results/paired_end/paired_end_R1_on_minus_antisense/ALFA/paired_end_R1_on_minus_antisense_Signal.UniqueMultiple.out.plus.bg
+8e23d52d7f635d927e292174f33168eb  results/paired_end/paired_end_R1_on_minus_sense/ALFA/paired_end_R1_on_minus_sense_Signal.UniqueMultiple.out.minus.bg
+5e90c760710980f4f4866dbe9aa32c6c  results/paired_end/paired_end_R1_on_minus_sense/ALFA/paired_end_R1_on_minus_sense_Signal.UniqueMultiple.out.plus.bg
+8e23d52d7f635d927e292174f33168eb  results/paired_end/paired_end_R1_on_plus_antisense/ALFA/paired_end_R1_on_plus_antisense_Signal.UniqueMultiple.out.minus.bg
+5e90c760710980f4f4866dbe9aa32c6c  results/paired_end/paired_end_R1_on_plus_antisense/ALFA/paired_end_R1_on_plus_antisense_Signal.UniqueMultiple.out.plus.bg
+5e90c760710980f4f4866dbe9aa32c6c  results/paired_end/paired_end_R1_on_plus_sense/ALFA/paired_end_R1_on_plus_sense_Signal.UniqueMultiple.out.minus.bg
+8e23d52d7f635d927e292174f33168eb  results/paired_end/paired_end_R1_on_plus_sense/ALFA/paired_end_R1_on_plus_sense_Signal.UniqueMultiple.out.plus.bg
+e5959524a2daf35da9249fb313920315  results/paired_end/paired_end_R1_on_minus_antisense/ALFA/paired_end_R1_on_minus_antisense.ALFA_feature_counts.tsv
+c406c1800d3690dd774aaa7e3c190523  results/paired_end/paired_end_R1_on_minus_sense/ALFA/paired_end_R1_on_minus_sense.ALFA_feature_counts.tsv
+c406c1800d3690dd774aaa7e3c190523  results/paired_end/paired_end_R1_on_plus_antisense/ALFA/paired_end_R1_on_plus_antisense.ALFA_feature_counts.tsv
+e5959524a2daf35da9249fb313920315  results/paired_end/paired_end_R1_on_plus_sense/ALFA/paired_end_R1_on_plus_sense.ALFA_feature_counts.tsv