Skip to content
Snippets Groups Projects
Commit bc168805 authored by Dominik Burri's avatar Dominik Burri
Browse files

Merge branch 'dev' into results_dir 

Update feature branch for later merging.
parents 9308bd60 f67aafb2
No related branches found
No related tags found
1 merge request!78Remove unnecessary files in results directory
Pipeline #11938 passed
Stage: test
Expand all Hide whitespace changes
Inline Side-by-side
Snakefile
+ 106
− 6
View file @ bc168805
......@@ -86,7 +86,6 @@ rule finish:
"summary_kallisto",
"genes_tpm.tsv")
rule start:
'''
Get samples
......@@ -808,10 +807,14 @@ rule kallisto_merge_genes:
gtf = get_sample('gtf')
output:
gn_out = os.path.join(
gn_tpm = os.path.join(
config["output_dir"],
"summary_kallisto",
"genes_tpm.tsv")
"genes_tpm.tsv"),
gn_counts = os.path.join(
config["output_dir"],
"summary_kallisto",
"genes_counts.tsv")
params:
dir_out = os.path.join(
......@@ -874,10 +877,14 @@ rule kallisto_merge_transcripts:
for i in pd.unique(samples_table.index.values)]),
output:
tx_out = os.path.join(
tx_tpm = os.path.join(
config["output_dir"],
"summary_kallisto",
"transcripts_tpm.tsv")
"transcripts_tpm.tsv"),
tx_counts = os.path.join(
config["output_dir"],
"summary_kallisto",
"transcripts_counts.tsv")
params:
dir_out = os.path.join(
......@@ -917,6 +924,87 @@ rule kallisto_merge_transcripts:
1> {log.stdout} 2> {log.stderr}"
rule pca_salmon:
input:
tpm = os.path.join(
config["output_dir"],
"summary_salmon",
"quantmerge",
"{molecule}_tpm.tsv"),
params:
tpm_filter = "0",
tpm_pseudocount = "1"
output:
out = directory(os.path.join(
config["output_dir"],
"zpca",
"pca_salmon_{molecule}"))
log:
stderr = os.path.join(
config["log_dir"],
"pca_salmon_{molecule}.stderr.log"),
stdout = os.path.join(
config["log_dir"],
"pca_salmon_{molecule}.stdout.log")
threads: 1
singularity:
"docker://zavolab/zpca:0.8"
shell:
"(zpca-tpm \
--tpm {input.tpm} \
--tpm-filter {params.tpm_filter} \
--tpm-pseudocount {params.tpm_pseudocount} \
--out {output.out} \
--verbose) \
1> {log.stdout} 2> {log.stderr}"
rule pca_kallisto:
input:
tpm = os.path.join(
config["output_dir"],
"summary_kallisto",
"{molecule}_tpm.tsv")
params:
tpm_filter = "0",
tpm_pseudocount = "1"
output:
out = directory(os.path.join(
config["output_dir"],
"zpca",
"pca_kallisto_{molecule}"))
log:
stderr = os.path.join(
config["log_dir"],
"pca_kallisto_{molecule}.stderr.log"),
stdout = os.path.join(
config["log_dir"],
"pca_kallisto_{molecule}.stdout.log")
threads: 1
singularity:
"docker://zavolab/zpca:0.8"
shell:
"(zpca-tpm \
--tpm {input.tpm} \
--tpm-filter {params.tpm_filter} \
--tpm-pseudocount {params.tpm_pseudocount} \
--out {output.out} \
--verbose) \
1> {log.stdout} 2> {log.stderr}"
rule star_rpm:
'''
Create stranded bedgraph coverage with STARs RPM normalisation
......@@ -1238,7 +1326,6 @@ rule alfa_qc:
-s {params.alfa_orientation}) &> {log}"
# cd {params.out_dir};
rule alfa_qc_all_samples:
'''
Run ALFA from stranded bedgraph files on all samples
......@@ -1405,6 +1492,19 @@ rule multiqc_report:
"ALFA",
"ALFA_plots_mqc.png"),
zpca_salmon = expand(os.path.join(
config["output_dir"],
"zpca",
"pca_salmon_{molecule}"),
molecule=["genes", "transcripts"]),
zpca_kallisto = expand(os.path.join(
config["output_dir"],
"zpca",
"pca_kallisto_{molecule}"),
molecule=["genes", "transcripts"]
),
multiqc_config = os.path.join(
config["output_dir"],
"multiqc_config.yaml")
......
images/dag_test_workflow.svg
+ 858
− 734
View file @ bc168805
This diff is collapsed.
images/rule_graph.svg
+ 439
− 391
View file @ bc168805
This diff is collapsed.
pipeline_documentation.md
+ 27
− 2
View file @ bc168805
......@@ -32,6 +32,8 @@ on installation and usage please see [here](README.md).
- [**salmon_quantmerge_transcripts**](#salmon_quantmerge_transcripts)
- [**kallisto_merge_genes**](#kallisto_merge_genes)
- [**kallisto_merge_transcripts**](#kallisto_merge_transcripts)
- [**pca_kallisto**](#pca_kallisto)
- [**pca_salmon**](#pca_salmon)
- [**generate_alfa_index**](#generate_alfa_index)
- [**alfa_qc**](#alfa_qc)
- [**alfa_qc_all_samples**](#alfa_qc_all_samples)
......@@ -400,7 +402,7 @@ Merge gene-level expression estimates for all samples with
- Gene TPM table (custom `.tsv`); used in
[**multiqc_report**](#multiqc_report)
- Gene read count table (custom `.tsv`); used in
[**multiqc_report**](#multiqc_report)
[**pca_salmon**](#pca_salmon) and [**pca_kallisto**](#pca_kallisto)
#### `salmon_quantmerge_transcripts`
......@@ -416,7 +418,7 @@ Merge transcript-level expression estimates for all samples with
- Transcript TPM table (custom `.tsv`); used in
[**multiqc_report**](#multiqc_report)
- Transcript read count table (custom `.tsv`); used in
[**multiqc_report**](#multiqc_report)
[**pca_salmon**](#pca_salmon) and [**pca_kallisto**](#pca_kallisto)
#### `kallisto_merge_genes`
......@@ -448,6 +450,28 @@ Merge transcript-level expression estimates for all samples with
- Transcript TPM table (custom `.tsv`)
- Transcript read count table (custom `.tsv`)
#### `pca_kallisto`
Run PCA analysis on kallisto genes and transcripts with [custom script][custom-script-zpca].
> Rule is run one time for transcript estimates and one time for genes estimates
- **Input**
- Transcript/Genes TPM table (custom `.tsv`)
- **Output**
- Directory with PCA plots, scree plot and top loading scores.
#### `pca_salmon`
Run PCA analysis on salmon genes and transcripts with [custom script][custom-script-zpca].
> Rule is run one time for transcript estimates and one time for genes estimates
- **Input**
- Transcript/Genes TPM table (custom `.tsv`)
- **Output**
- Directory with PCA plots, scree plot and top loading scores.
#### `generate_alfa_index`
Create index for [**ALFA**](#third-party-software-used).
......@@ -714,6 +738,7 @@ Generate pseudoalignments of reads to transcripts with
[custom-script-gtf-to-bed12]: <https://git.scicore.unibas.ch/zavolan_group/tools/gtf_transcript_type_to_bed12>
[custom-script-tin]: <https://git.scicore.unibas.ch/zavolan_group/tools/tin_score_calculation>
[custom-script-merge-kallisto]: <https://github.com/zavolanlab/merge_kallisto>
[custom-script-zpca]: <https://github.com/zavolanlab/zpca>
[docs-alfa]: <https://github.com/biocompibens/ALFA#manual>
[docs-bedgraphtobigwig]: <http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/>
[docs-bedtools]: <https://bedtools.readthedocs.io/en/latest/>
......
scripts/fg_extract_transcripts.py deleted 100755 → 0
+ 0
− 104
View file @ 9308bd60
# _____________________________________________________________________________
# -----------------------------------------------------------------------------
# import needed (external) modules
# -----------------------------------------------------------------------------
try:
import HTSeq
except(Exception):
raise("[ERROR] HTSeq was not imported properly. Exiting.")
sys.exit(-1)
try:
from argparse import ArgumentParser, RawTextHelpFormatter, FileType
except(Exception):
raise("[ERROR] argparse was not imported properly. Exiting.")
sys.exit(-1)
try:
import sys
except(Exception):
raise("[ERROR] sys was not imported properly. Exiting.")
sys.exit(-1)
# _____________________________________________________________________________
# -----------------------------------------------------------------------------
# Main function
# -----------------------------------------------------------------------------
def main():
""" Main function """
__doc__ = "Filter gtf file based on the transcript biotype."
parser = ArgumentParser(description=__doc__, formatter_class=RawTextHelpFormatter)
parser.add_argument("--gtf",
dest = "gtf",
help = "Annotation file in GTF format with transcript support level information",
required = True,
metavar = "FILE")
parser.add_argument("--out",
dest = "out",
help = "GTF output file",
required = True,
metavar="FILE")
parser.add_argument("--transcript_biotype",
dest = "transcript_biotype",
help = "Transcript biotype or transcript type to include",
required = True,
default = "protein_coding,lincRNA")
parser.add_argument("-v",
"--verbose",
action = "store_true",
dest = "verbose",
default = False,
required = False,
help = "Verbose")
# _________________________________________________________________________
# -------------------------------------------------------------------------
# get the arguments
# -------------------------------------------------------------------------
try:
options = parser.parse_args()
except(Exception):
parser.print_help()
if len(sys.argv) == 1:
parser.print_help()
sys.exit(1)
if options.verbose:
sys.stdout.write("Parsing gtf file and keep %s transcripts\n" % (str(options.transcript_biotype)))
list_of_transcript_types = str(options.transcript_biotype).split(",")
# parse gtf file
gtf = HTSeq.GFF_Reader(options.gtf)
# open output file
w = open(options.out, 'w')
for gtf_line in gtf:
if "transcript_biotype" in gtf_line.attr:
if gtf_line.attr['transcript_biotype'] in list_of_transcript_types:
w.write(gtf_line.get_gff_line())
if "transcript_type" in gtf_line.attr:
if gtf_line.attr['transcript_type'] in list_of_transcript_types:
w.write(gtf_line.get_gff_line())
w.close()
# _____________________________________________________________________________
# -----------------------------------------------------------------------------
# Call the Main function and catch Keyboard interrups
# -----------------------------------------------------------------------------
if __name__ == '__main__':
try:
main()
except KeyboardInterrupt:
sys.stderr.write("User interrupt!\n")
sys.exit(0)
scripts/heatmap_and_clustermap.py deleted 100755 → 0
+ 0
− 112
View file @ 9308bd60
#!/usr/bin/env python
# _____________________________________________________________________________
# -----------------------------------------------------------------------------
# import needed (external) modules
# -----------------------------------------------------------------------------
import sys
import os
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from argparse import ArgumentParser, RawTextHelpFormatter
# _____________________________________________________________________________
# -----------------------------------------------------------------------------
# Main function
# -----------------------------------------------------------------------------
def main():
""" Main function """
__doc__ = "Create heatmap and clustermap based on a table (rows are genes/transcripts, columns are samples)."
parser = ArgumentParser(
description=__doc__,
formatter_class=RawTextHelpFormatter
)
parser.add_argument(
"--tpm",
dest="tpm",
help="TPM table (tsv)",
required=True,
metavar="FILE"
)
parser.add_argument(
"--out",
dest="out",
help="Output directory",
required=True,
metavar="FILE"
)
parser.add_argument(
"-v",
"--verbose",
action="store_true",
dest="verbose",
default=False,
required=False,
help="Verbose"
)
# _________________________________________________________________________
# -------------------------------------------------------------------------
# get the arguments
# -------------------------------------------------------------------------
try:
options = parser.parse_args()
except(Exception):
parser.print_help()
if len(sys.argv) == 1:
parser.print_help()
sys.exit(1)
if options.verbose:
sys.stdout.write(f"Creating output directory: {options.out} {os.linesep}")
if not os.path.exists(options.out):
os.makedirs(options.out)
if options.verbose:
sys.stdout.write(f"Reading: {options.tpm} {os.linesep}")
df = pd.read_csv(options.tpm, header=0, sep="\t")
df.set_index(["Name"], inplace=True)
df = df + 1/32
df = np.log2(df)
if options.verbose:
sys.stdout.write(f"Generating heatmap {os.linesep}")
plt.rcParams["figure.figsize"] = (20,20)
g = sns.heatmap(df.corr(method="pearson"))
plt.savefig(os.path.join(options.out, "heatmap.pdf"))
if options.verbose:
sys.stdout.write(f"Generating clustermap {os.linesep}")
plt.rcParams["figure.figsize"] = (20,20)
g = sns.clustermap(df.corr(method="pearson"))
plt.savefig(os.path.join(options.out, "clustermap.pdf"))
if options.verbose:
sys.stdout.write(f"Done {os.linesep}")
# _____________________________________________________________________________
# -----------------------------------------------------------------------------
# Call the Main function and catch Keyboard interrups
# -----------------------------------------------------------------------------
if __name__ == '__main__':
try:
main()
except KeyboardInterrupt:
sys.stderr.write("User interrupt!" + os.linesep)
sys.exit(0)
scripts/perform_PCA.py deleted 100755 → 0
+ 0
− 210
View file @ 9308bd60
#!/usr/bin/env python
# _____________________________________________________________________________
# -----------------------------------------------------------------------------
# import needed (external) modules
# -----------------------------------------------------------------------------
import sys
import os
import pandas as pd
import numpy as np
import random as rd
from sklearn.decomposition import PCA
from sklearn import preprocessing
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from argparse import ArgumentParser, RawTextHelpFormatter
# _____________________________________________________________________________
# -----------------------------------------------------------------------------
# Main function
# -----------------------------------------------------------------------------
def main():
""" Main function """
__doc__ = "Perform PCA based on a table (rows are genes/transcripts, columns are samples)."
parser = ArgumentParser(
description=__doc__,
formatter_class=RawTextHelpFormatter
)
parser.add_argument(
"--tpm",
dest="tpm",
help="TPM table (tsv)",
required=True,
metavar="FILE"
)
parser.add_argument(
"--out",
dest="out",
help="Output directory",
required=True,
metavar="FILE"
)
parser.add_argument(
"-v",
"--verbose",
action="store_true",
dest="verbose",
default=False,
required=False,
help="Verbose"
)
# _________________________________________________________________________
# -------------------------------------------------------------------------
# get the arguments
# -------------------------------------------------------------------------
try:
options = parser.parse_args()
except(Exception):
parser.print_help()
if len(sys.argv) == 1:
parser.print_help()
sys.exit(1)
if options.verbose:
sys.stdout.write(f"Creating output directory: {options.out} {os.linesep}")
if not os.path.exists(options.out):
os.makedirs(options.out)
if options.verbose:
sys.stdout.write(f"Reading: {options.tpm} {os.linesep}")
df = pd.read_csv(options.tpm, header=0, sep="\t")
df.set_index(["Name"], inplace=True)
if options.verbose:
sys.stdout.write(f"Performing PCA {os.linesep}")
scaled_data = preprocessing.scale(df.T)
pca = PCA()
pca.fit(scaled_data)
pca_data = pca.transform(scaled_data)
per_var = np.round(pca.explained_variance_ratio_* 100, decimals=1)
labels = ['PC' + str(x) for x in range(1, len(per_var)+1)]
plt.figure()
plt.rcParams["figure.figsize"] = (20,20)
plt.bar(x=range(1,len(per_var)+1), height=per_var, tick_label=labels)
plt.ylabel('Percentage of Explained Variance')
plt.xlabel('Principal Component')
plt.title('Scree Plot')
plt.xticks(rotation='vertical')
plt.savefig(os.path.join(options.out, "scree_plot.pdf"))
plt.close()
pca_df = pd.DataFrame(pca_data, index=[*df], columns=labels)
if options.verbose:
sys.stdout.write(f"Generating plots in: {options.out} {os.linesep}")
# -------------------------------------------------------------------------
# PCA 1st and 2nd component
# -------------------------------------------------------------------------
plt.figure()
plt.rcParams["figure.figsize"] = (20,20)
plt.scatter(pca_df.PC1, pca_df.PC2)
plt.xlabel('PC1 - {0}%'.format(per_var[0]), fontsize=22)
plt.ylabel('PC2 - {0}%'.format(per_var[1]), fontsize=22)
plt.title('PCA components 1 & 2', fontsize=22)
for sample in pca_df.index:
plt.annotate(sample, (pca_df.PC1.loc[sample], pca_df.PC2.loc[sample]))
plt.savefig(os.path.join(options.out, "PCA_1_2.pdf"))
plt.close()
# -------------------------------------------------------------------------
# PCA 2nd and 3rd component
# -------------------------------------------------------------------------
plt.figure()
plt.rcParams["figure.figsize"] = (20,20)
plt.scatter(pca_df.PC2, pca_df.PC3)
plt.xlabel('PC2 - {0}%'.format(per_var[1]), fontsize=22)
plt.ylabel('PC3 - {0}%'.format(per_var[2]), fontsize=22)
plt.title('PCA components 2 & 3', fontsize=22)
for sample in pca_df.index:
plt.annotate(sample, (pca_df.PC2.loc[sample], pca_df.PC3.loc[sample]))
plt.savefig(os.path.join(options.out, "PCA_2_3.pdf"))
plt.close()
# -------------------------------------------------------------------------
# PCA 1st and 3rd component
# -------------------------------------------------------------------------
plt.figure()
plt.rcParams["figure.figsize"] = (20,20)
plt.scatter(pca_df.PC1, pca_df.PC3)
plt.xlabel('PC1 - {0}%'.format(per_var[0]), fontsize=22)
plt.ylabel('PC3 - {0}%'.format(per_var[2]), fontsize=22)
plt.title('PCA components 1 & 3', fontsize=22)
for sample in pca_df.index:
plt.annotate(sample, (pca_df.PC1.loc[sample], pca_df.PC3.loc[sample]))
plt.savefig(os.path.join(options.out, "PCA_1_3.pdf"))
plt.close()
# -------------------------------------------------------------------------
# PCA 3D
# -------------------------------------------------------------------------
labels = []
pc1 = []
pc2 = []
pc3 = []
for i, row in pca_df.iterrows():
labels.append(i)
pc1.append(row["PC1"])
pc2.append(row["PC2"])
pc3.append(row["PC3"])
fig = plt.figure(figsize=(20,20))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(pc1, pc2, pc3)
ax.set_xlabel('PC1 - {0}%'.format(per_var[0]))
ax.set_ylabel('PC2 - {0}%'.format(per_var[1]))
ax.set_zlabel('PC3 - {0}%'.format(per_var[2]))
for label, x, y, z in zip(labels, pc1, pc2, pc3):
ax.text(x, y, z, label)
plt.savefig(os.path.join(options.out, "PCA_3D.pdf"))
# -------------------------------------------------------------------------
# Write loading scores
# -------------------------------------------------------------------------
if options.verbose:
sys.stdout.write(f"Writing loading scores in: {options.out} {os.linesep}")
loading_scores = pd.Series(pca.components_[0], index=list(df.index))
sorted_loading_scores = loading_scores.abs().sort_values(ascending=False)
sorted_loading_scores.to_csv(os.path.join(options.out, "loading_scores.tsv"), header=False, sep="\t")
if options.verbose:
sys.stdout.write(f"Done {os.linesep}")
# _____________________________________________________________________________
# -----------------------------------------------------------------------------
# Call the Main function and catch Keyboard interrups
# -----------------------------------------------------------------------------
if __name__ == '__main__':
try:
main()
except KeyboardInterrupt:
sys.stderr.write("User interrupt!" + os.linesep)
sys.exit(0)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment