diff --git a/Snakefile b/Snakefile
index 6925a076dd3f0af191d7ed42128c6f436aef480d..c98e58ec862c32856ee4ed22e03cf0550377f378 100644
--- a/Snakefile
+++ b/Snakefile
@@ -1,6 +1,6 @@
#---------------------------
# @author: Kedi Cao
-# @email: ckedi.cao@ukr.de
+# @email: kedi.cao@ukr.de
# @date: 28.08.2024
#---------------------------
@@ -36,12 +36,7 @@ kraken_lib = config["kraken_lib"]
# REMAP = config["REMAP"]
-import os
-def get_results_filename(filenames):
- paths = [os.path.join(RESULTS_DIR, filename) for filename in filenames]
- existing_paths = [path if os.path.exists(path) else [] for path in paths]
- return existing_paths
################################################################################
input_files = [
@@ -60,10 +55,18 @@ if config["remove_contamination"]=="no":
elif config["remove_contamination"]=="yes":
input_files.append(expand("{p}/{s}/{f}",p=WORK_DIR, s=SAMPLES, f=["single_clean.fastq","paired_1_clean.fastq","paired_2_clean.fastq"]))
input_files.append(expand("{p}/{s}/remapped_non_ref.bam", p=WORK_DIR, s=SAMPLES))
- # input_files.append(RESULTS_DIR + "/read_numbers.png")
if config["ANALYSIS"]=="yes":
- input_files.append(expand("{p}/{f}", p=RESULTS_DIR, f=["read_numbers.png","analysis_table.png","gc_content_dist.png","contig_length_dist.png","heatmap_coverage.png"]))
+ input_files.append(expand("{p}/{f}", p=RESULTS_DIR, f=["analysis_table.png","gc_content_dist.png","contig_length_dist.png","heatmap_coverage.png"]))
+ if config["SICKLE"]=="no":
+ input_files.append(expand("{p}/{f}", p=RESULTS_DIR, f=["read_numbers.png"]))
+ elif config["SICKLE"]=="yes":
+ if config["remove_contamination"]=="no":
+ input_files.append(expand("{p}/{f}", p=RESULTS_DIR, f=["read_numbers_before_filter.png","read_numbers_after_filter.png"]))
+ elif config["remove_contamination"]=="yes":
+ input_files.append(expand("{p}/{f}", p=RESULTS_DIR, f=["read_numbers_before_clean_filter.png","read_numbers_after_clean.png","read_numbers_after_clean_filter.png"]))
+
+
################################################################################
diff --git a/modules/analysis.smk b/modules/analysis.smk
index 404842875517f9107994a94e620a0ec20f5ef8dc..5204b0b237504272322ea0a379717c5862e2bd2b 100644
--- a/modules/analysis.smk
+++ b/modules/analysis.smk
@@ -1,21 +1,107 @@
-rule unmapped_analysis:
- input:
- single = expand("{p}/{s}/single.fastq", p=WORK_DIR, s=SAMPLES),
- paired1 = expand("{p}/{s}/paired.1.fastq", p=WORK_DIR, s=SAMPLES),
- paired2 = expand("{p}/{s}/paired.2.fastq", p=WORK_DIR, s=SAMPLES)
- output:
- RESULTS_DIR + "/read_numbers.png"
- conda:
- "envs/notebooks.yml"
- log:
- "logs/analysis/unmapped_analysis.log"
- benchmark:
- "benchmarks/analysis/unmapped_analysis.txt"
- notebook:
- MYPATH + "/modules/notebooks/read_numbers_in_bar_graph.py.ipynb"
+if config["SICKLE"]=="no":
+ rule unmapped_analysis:
+ input:
+ single = expand("{p}/{s}/single.fastq", p=WORK_DIR, s=SAMPLES),
+ paired1 = expand("{p}/{s}/paired.1.fastq", p=WORK_DIR, s=SAMPLES),
+ paired2 = expand("{p}/{s}/paired.2.fastq", p=WORK_DIR, s=SAMPLES)
+ output:
+ RESULTS_DIR + "/read_numbers.png"
+ conda:
+ "envs/notebooks.yml"
+ log:
+ "logs/analysis/unmapped_analysis.log"
+ benchmark:
+ "benchmarks/analysis/unmapped_analysis.txt"
+ notebook:
+ MYPATH + "/modules/notebooks/read_numbers_in_bar_graph.py.ipynb"
+
+elif config["SICKLE"]=="yes":
+ if config["remove_contamination"]=="no":
+ rule unmapped_analysis_before_filter:
+ input:
+ single = expand("{p}/{s}/single.fastq", p=WORK_DIR, s=SAMPLES),
+ paired1 = expand("{p}/{s}/paired.1.fastq", p=WORK_DIR, s=SAMPLES),
+ paired2 = expand("{p}/{s}/paired.2.fastq", p=WORK_DIR, s=SAMPLES)
+ output:
+ RESULTS_DIR + "/read_numbers_before_filter.png"
+ conda:
+ "envs/notebooks.yml"
+ log:
+ "logs/analysis/unmapped_analysis_before_filter.log"
+ benchmark:
+ "benchmarks/analysis/unmapped_analysis_before_filter.txt"
+ notebook:
+ MYPATH + "/modules/notebooks/read_numbers_in_bar_graph.py.ipynb"
+
+ rule unmapped_analysis_after_filter:
+ input:
+ single = expand("{p}/{s}/sickle.single.fastq", p=WORK_DIR, s=SAMPLES),
+ paired1 = expand("{p}/{s}/sickle.paired.1.fastq", p=WORK_DIR, s=SAMPLES),
+ paired2 = expand("{p}/{s}/sickle.paired.2.fastq", p=WORK_DIR, s=SAMPLES)
+ output:
+ RESULTS_DIR + "/read_numbers_after_filter.png"
+ conda:
+ "envs/notebooks.yml"
+ log:
+ "logs/analysis/unmapped_analysis_after_filter.log"
+ benchmark:
+ "benchmarks/analysis/unmapped_analysis_after_filter.txt"
+ notebook:
+ MYPATH + "/modules/notebooks/read_numbers_in_bar_graph.py.ipynb"
+
+ elif config["remove_contamination"]=="yes":
+ rule unmapped_analysis_before_clean_filter:
+ input:
+ single = expand("{p}/{s}/single.fastq", p=WORK_DIR, s=SAMPLES),
+ paired1 = expand("{p}/{s}/paired.1.fastq", p=WORK_DIR, s=SAMPLES),
+ paired2 = expand("{p}/{s}/paired.2.fastq", p=WORK_DIR, s=SAMPLES)
+ output:
+ RESULTS_DIR + "/read_numbers_before_clean_filter.png"
+ conda:
+ "envs/notebooks.yml"
+ log:
+ "logs/analysis/unmapped_analysis_before_clean_filter.log"
+ benchmark:
+ "benchmarks/analysis/unmapped_analysis_before_clean_filte.txt"
+ notebook:
+ MYPATH + "/modules/notebooks/read_numbers_in_bar_graph.py.ipynb"
+
+ rule unmapped_analysis_after_clean:
+ input:
+ single = expand("{p}/{s}/single_clean.fastq", p=WORK_DIR, s=SAMPLES),
+ paired1 = expand("{p}/{s}/paired_1_clean.fastq", p=WORK_DIR, s=SAMPLES),
+ paired2 = expand("{p}/{s}/paired_2_clean.fastq", p=WORK_DIR, s=SAMPLES)
+ output:
+ RESULTS_DIR + "/read_numbers_after_clean.png"
+ conda:
+ "envs/notebooks.yml"
+ log:
+ "logs/analysis/unmapped_analysis_after_clean.log"
+ benchmark:
+ "benchmarks/analysis/unmapped_analysis_after_cleaning.txt"
+ notebook:
+ MYPATH + "/modules/notebooks/read_numbers_in_bar_graph.py.ipynb"
+
+ rule unmapped_analysis_after_clean_filter:
+ input:
+ single = expand("{p}/{s}/sickle.single_clean.fastq", p=WORK_DIR, s=SAMPLES),
+ paired1 = expand("{p}/{s}/sickle.paired_1_clean.fastq", p=WORK_DIR, s=SAMPLES),
+ paired2 = expand("{p}/{s}/sickle.paired_2_clean.fastq", p=WORK_DIR, s=SAMPLES)
+ output:
+ RESULTS_DIR + "/read_numbers_after_clean_filter.png"
+ conda:
+ "envs/notebooks.yml"
+ log:
+ "logs/analysis/unmapped_analysis_after_clean_filter.log"
+ benchmark:
+ "benchmarks/analysis/unmapped_analysis_after_clean_filter.txt"
+ notebook:
+ MYPATH + "/modules/notebooks/read_numbers_in_bar_graph.py.ipynb"
+
+
+
-
rule assembly_analysis_table:
input:
contigs = expand("{p}/{s}/{assemblr}.contigs.fa", p=WORK_DIR, s=SAMPLES, assemblr=ASSEMBLER)
@@ -46,6 +132,22 @@ rule assembly_analysis:
notebook:
MYPATH + "/modules/notebooks/content_distribution.py.ipynb"
+
+rule calculate_coverage:
+ input:
+ bams = WORK_DIR + "/{sample}/non_ref_new.bam",
+ bais = WORK_DIR + "/{sample}/non_ref_new.bam.bai"
+ output:
+ cov = WORK_DIR + "/{sample}/coverage.txt"
+ conda:
+ "envs/samtools.yml"
+ log:
+ "logs/analysis/{sample}_calculate_coverage.log"
+ benchmark:
+ "benchmarks/analysis/{sample}_calculate_coverage.txt"
+ shell:
+ "samtools coverage {input.bams} > {output.cov}"
+
rule coverage_analysis:
input:
@@ -61,20 +163,3 @@ rule coverage_analysis:
"benchmarks/analysis/coverage_analysis.txt"
notebook:
MYPATH + "/modules/notebooks/contig_coverage_heatmap.py.ipynb"
-
-# rule coverage_analysis:
-# input:
-# cov = RESULTS_DIR + "/{sample}_coverage.txt"
-# bams = expand(WORK_DIR + "/{s}/non_ref_new.bam", s=SAMPLES),
-# bais = expand(WORK_DIR + "/{s}/non_ref_new.bam.bai", s=SAMPLES),
-# fasta = RESULTS_DIR + "/supercontigs.fa"
-# output:
-# RESULTS_DIR + "/heatmap_coverage.png"
-# conda:
-# "envs/notebooks.yml"
-# log:
-# "logs/analysis/coverage_analysis.log"
-# benchmark:
-# "benchmarks/analysis/coverage_analysis.txt"
-# notebook:
-# MYPATH + "/modules/notebooks/contig_coverage_heatmap.py.ipynb"
\ No newline at end of file
diff --git a/modules/contigmap.smk b/modules/contigmap.smk
index 4e6a46e2427a1023841f9a4b05df442b7b5ea8b7..9b89be357d442843c7872441398813b8cf1df66c 100644
--- a/modules/contigmap.smk
+++ b/modules/contigmap.smk
@@ -165,20 +165,4 @@ rule index_sorted:
benchmark:
"benchmarks/contigmap/{sample}_index_sorted.txt"
shell:
- "samtools index {input}"
-
-
-rule calculate_coverage:
- input:
- bams = WORK_DIR + "/{sample}/non_ref_new.bam",
- bais = WORK_DIR + "/{sample}/non_ref_new.bam.bai"
- output:
- cov = WORK_DIR + "/{sample}/coverage.txt"
- conda:
- "envs/samtools.yml"
- log:
- "logs/contigmap/{sample}_calculate_coverage.log"
- benchmark:
- "benchmarks/contigmap/{sample}_calculate_coverage.txt"
- shell:
- "samtools coverage {input.bams} > {output.cov}"
\ No newline at end of file
+ "samtools index {input}"
\ No newline at end of file
diff --git a/modules/crop_remapped.smk b/modules/crop_remapped.smk
index 8690fe671c717d30196e35ff505d9d7ed9b8bf6d..bea5f1dea227632a678852cfd53d0a84b3894116 100644
--- a/modules/crop_remapped.smk
+++ b/modules/crop_remapped.smk
@@ -70,7 +70,7 @@ elif config["SICKLE"]=="yes":
rule popins2_sort:
input:
- os.path.join(RESULTS_DIR, "read_numbers.png") if config["ANALYSIS"]=="yes" else [],
+ # os.path.join(RESULTS_DIR, "read_numbers.png") if config["ANALYSIS"]=="yes" else [],
mates = WORK_DIR + "/{sample}/mates.bam"
output:
WORK_DIR + "/{sample}/non_ref.bam"
@@ -89,7 +89,7 @@ elif config["SICKLE"]=="yes":
rule popins2_sickle:
input:
- os.path.join(RESULTS_DIR, "read_numbers.png") if config["ANALYSIS"]=="yes" else [],
+ # os.path.join(RESULTS_DIR, "read_numbers.png") if config["ANALYSIS"]=="yes" else [],
single = WORK_DIR + "/{sample}/single_clean.fastq",
pair_1 = WORK_DIR + "/{sample}/paired_1_clean.fastq",
pair_2 = WORK_DIR + "/{sample}/paired_2_clean.fastq"
diff --git a/modules/crop_unmapped.smk b/modules/crop_unmapped.smk
index 6e95d5391af0486b3c41b9d23a33ca4f1ea8d03b..036c5c37f13087f6c3887ce7d0a89f0541512ab1 100644
--- a/modules/crop_unmapped.smk
+++ b/modules/crop_unmapped.smk
@@ -68,7 +68,7 @@ elif config["SICKLE"]=="yes":
rule popins2_sort:
input:
- os.path.join(RESULTS_DIR, "read_numbers.png") if config["ANALYSIS"]=="yes" else [],
+ # os.path.join(RESULTS_DIR, "read_numbers.png") if config["ANALYSIS"]=="yes" else [],
mates = WORK_DIR + "/{sample}/mates.bam"
output:
WORK_DIR + "/{sample}/non_ref.bam"
@@ -87,7 +87,7 @@ elif config["SICKLE"]=="yes":
rule popins2_sickle:
input:
- os.path.join(RESULTS_DIR, "read_numbers.png") if config["ANALYSIS"]=="yes" else [],
+ # os.path.join(RESULTS_DIR, "read_numbers.png") if config["ANALYSIS"]=="yes" else [],
single = WORK_DIR + "/{sample}/single.fastq",
pair_1 = WORK_DIR + "/{sample}/paired.1.fastq",
pair_2 = WORK_DIR + "/{sample}/paired.2.fastq"
diff --git a/modules/genotype.smk b/modules/genotype.smk
index f21caa67a6a1f54ed3f1b9a994b579ab41c434c2..df1c581ca1b73f8252471e9a103e9234c683b816 100644
--- a/modules/genotype.smk
+++ b/modules/genotype.smk
@@ -1,11 +1,7 @@
-# RESULTS_DIR = config["RESULTS_DIR"]
-# WORK_DIR = config["WORK_DIR"]
-# SAMPLES = config["SAMPLES"]
-
rule sort_vcf:
input:
RESULTS_DIR + "/place-finish.done",
- ins=rules.popins2_place_refalign.output.vcf
+ ins=RESULTS_DIR + "/insertions_unsorted.vcf"
output:
temp(RESULTS_DIR + "/insertions.vcf")
conda:
@@ -15,7 +11,6 @@ rule sort_vcf:
benchmark:
"benchmarks/genotype/sort_vcf.txt"
shell:
- "bcftools view {input.ins} | "
"bcftools sort {input.ins} "
" --output-file {output}"
" --output-type v"
@@ -33,7 +28,7 @@ rule popins2_genotype:
shell:
"{POPINS2_BIN} genotype "
" --prefix {WORK_DIR} "
- " --insertions {rules.sort_vcf.output} "
+ " --insertions {input} "
" --contigs {rules.popins2_merge_contigs.output.supercontigs} "
" --reference {REFERENCE} "
" {wildcards.sample} "
diff --git a/modules/notebooks/contig_coverage_heatmap.py.ipynb b/modules/notebooks/contig_coverage_heatmap.py.ipynb
index 20ae7d97088991d7b4637e0cbd053d4475c74899..46ac0ca7995116e6cab2739d39d7648b949166ba 100644
--- a/modules/notebooks/contig_coverage_heatmap.py.ipynb
+++ b/modules/notebooks/contig_coverage_heatmap.py.ipynb
@@ -126,20 +126,8 @@
" 'contig_id':contig_id,\n",
" 'contig_length':contig_length\n",
"})\n",
- "df_sort_contig\n",
"df_sort_contig = df_sort_contig.sort_values(by=['contig_length'], ascending=True)\n",
- "order_by_length = df_sort_contig['contig_id'].tolist()\n",
- "order_by_length"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 89,
- "id": "39150e4b-d2f4-43d4-a6da-6427453899fa",
- "metadata": {},
- "outputs": [],
- "source": [
- "# df = pd.read_csv(coverage_files[0], sep='\\t')"
+ "order_by_length = df_sort_contig['contig_id'].tolist()"
]
},
{
@@ -167,12 +155,7 @@
" coverage_df = df\n",
" else:\n",
" # Merge with the main DataFrame on the contig name\n",
- " coverage_df = pd.merge(coverage_df, df, on='#rname', how='outer')\n",
- " # Set the contig names as the index\n",
- "# coverage_df.set_index('#rname', inplace=True)\n",
- "# Replace NaN values with zeros (if any contigs are missing in some samples)\n",
- "# coverage_df.fillna(0, inplace=True)\n",
- "# coverage_df.fillna(0)"
+ " coverage_df = pd.merge(coverage_df, df, on='#rname', how='outer')\n"
]
},
{
@@ -187,6 +170,16 @@
"coverage_df = coverage_df.reindex(order_by_length)"
]
},
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ab7e6633",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "coverage_df_cleaned = calculator.remove_outliers_replace_with_bounds(coverage_df)"
+ ]
+ },
{
"cell_type": "code",
"execution_count": 85,
@@ -195,9 +188,8 @@
"outputs": [],
"source": [
"# Plot the heatmap\n",
- "plt.figure(figsize=(10, 10))\n",
- "sns.set(font_scale = 1.3)\n",
- "sns.heatmap(coverage_df, annot=False, cmap='YlOrBr')\n",
+ "sns.set(font_scale = 0.8)\n",
+ "sns.heatmap(coverage_df_cleaned, annot=False, cmap='YlOrBr')\n",
"plt.title('Coverage of Contigs')\n",
"plt.xlabel('Samples')\n",
"plt.ylabel('Contigs')\n",
@@ -206,19 +198,6 @@
"plt.savefig(snakemake.output[0], dpi=400, bbox_inches = 'tight')"
]
},
- {
- "cell_type": "code",
- "execution_count": 79,
- "id": "f0e2f837-0307-4269-928d-417d89303aeb",
- "metadata": {},
- "outputs": [],
- "source": [
- "\n",
- "\n",
- "\n",
- "# plt.savefig(wkdir+'/contig_coverage_no_outlier', dpi=400, bbox_inches = 'tight')"
- ]
- },
{
"cell_type": "code",
"execution_count": null,
diff --git a/modules/notebooks/src/stat_func/calculator.py b/modules/notebooks/src/stat_func/calculator.py
index 8c080fcc81bd3156607f33dcf7c74a093138806d..f85adca87a0f37575728d0d14620599bb0f6ee50 100644
--- a/modules/notebooks/src/stat_func/calculator.py
+++ b/modules/notebooks/src/stat_func/calculator.py
@@ -250,6 +250,39 @@ def replace_column_outliers_with_max(df, column):
return df
+def remove_outliers_replace_with_bounds(df):
+ """
+ Removes outliers from the 1st and 4th quartiles and replaces them with the lower and upper bounds.
+
+ Parameters:
+ df: DataFrame where the first column is the index and the next columns are the contig coverages for samples.
+
+ Returns:
+ df: DataFrame with outliers replaced by the respective bounds.
+ """
+
+ # Apply the process to each column starting from the second (ignoring index column)
+ for col in df.columns:
+ # Calculate Q1 (25th percentile) and Q3 (75th percentile)
+ Q1 = df[col].quantile(0.25)
+ Q3 = df[col].quantile(0.75)
+
+ # Calculate IQR
+ IQR = Q3 - Q1
+
+ # Define bounds for outliers
+ lower_bound = Q1 - 1.5 * IQR
+ upper_bound = Q3 + 1.5 * IQR
+
+ # Replace values below the lower bound with the lower bound
+ df[col] = df[col].apply(lambda x: lower_bound if x < lower_bound else x)
+
+ # Replace values above the upper bound with the upper bound
+ df[col] = df[col].apply(lambda x: upper_bound if x > upper_bound else x)
+
+ return df
+
+
#quick function for getting an average of a list of numbers as a string
def quickstats(x):
x_array = np.asarray(x,dtype=np.float64)
diff --git a/modules/position.smk b/modules/position.smk
index 6fcfed96c33ed85bab45426192e10e0c7dd53ea6..ae9cf8be2ef4ff65c38fa2424636ef6499823fed 100644
--- a/modules/position.smk
+++ b/modules/position.smk
@@ -1,10 +1,3 @@
-# WORK_DIR = config["WORK_DIR"]
-# RESULTS_DIR = config["RESULTS_DIR"]
-# REFERENCE = config["REFERENCE"]
-# SAMPLES = config["SAMPLES"]
-# ASSEMBLER = config["ASSEMBLER"]
-
-
if config["remove_contamination"] == 'yes':
@@ -54,7 +47,6 @@ rule popins2_merge_locations:
input:
os.path.join(RESULTS_DIR, "heatmap_coverage.png") if config["ANALYSIS"]=="yes" else [],
expand(WORK_DIR + "/{s}/locations.txt", s=SAMPLES)
- # RESULTS_DIR + "/heatmap_coverage.png"
output:
loc=RESULTS_DIR + "/locations.txt"
log: