From 8d5ddcb1bb2135532817f465966a8d95189bf386 Mon Sep 17 00:00:00 2001
From: mariabernard <maria.bernard@jouy.inra.fr>
Date: Wed, 19 Jun 2019 17:20:01 +0200
Subject: [PATCH] 1000RNASEQ chicken ASE: change bam index output name, use
.bam.bai for phASER and .bai for ASEReadCounter, change phASER output csv
file
---
Snakemake/1000RNASeq_chicken/ASE/Snakefile | 4 ++--
.../1000RNASeq_chicken/ASE/rules/ASE_counter.smk | 11 ++++++-----
Snakemake/1000RNASeq_chicken/ASE/rules/basics.smk | 2 +-
3 files changed, 9 insertions(+), 8 deletions(-)
diff --git a/Snakemake/1000RNASeq_chicken/ASE/Snakefile b/Snakemake/1000RNASeq_chicken/ASE/Snakefile
index be6201e..5a67cbf 100644
--- a/Snakemake/1000RNASeq_chicken/ASE/Snakefile
+++ b/Snakemake/1000RNASeq_chicken/ASE/Snakefile
@@ -138,14 +138,14 @@ for sample in table["sample_name"].unique():
# # counter
for sample in table["sample_name"].unique():
final_outputs.append('Results/Counter/' + sample + '_ASEReadCounter.csv')
- final_outputs.append('Results/Counter/' + sample + '_phASER.csv')
+ final_outputs.append('Results/Counter/' + sample + '_phASER.allelic_counts.txt')
# print(final_outputs)
# ['Results/genomeMasked/reference_masked.fa',
# 'Results/STAR_Aln_2/samplePE_rg_genomic_pp_uniq_rmdup_split.bam',
# 'Results/STAR_Aln_2/sampleSE_rg_genomic_uniq_rmdup_split.bam',
# 'Results/Counter/samplePE_ASEReadCounter.csv', 'Results/Counter/sampleSE_ASEReadCounter.csv',
-# 'Results/Counter/samplePE_phASER.csv', 'Results/Counter/sampleSE_phASER.csv']
+# 'Results/Counter/samplePE_phASER.allelic_counts.txt', 'Results/Counter/sampleSE_phASER.allelic_counts.txt']
rule all:
input:
diff --git a/Snakemake/1000RNASeq_chicken/ASE/rules/ASE_counter.smk b/Snakemake/1000RNASeq_chicken/ASE/rules/ASE_counter.smk
index 95cfa1f..e7ceb3d 100644
--- a/Snakemake/1000RNASeq_chicken/ASE/rules/ASE_counter.smk
+++ b/Snakemake/1000RNASeq_chicken/ASE/rules/ASE_counter.smk
@@ -28,10 +28,10 @@ def get_splitNCigar_bai(wildcards):
# paired end, filter on properly paired and rmdup
if len(sample_table["forward_read"]) == len(sample_table["reverse_read"].dropna()) :
- return "Results/STAR_Aln_2/" + wildcards.sample + "_rg_genomic_pp_rmdup_uniq_split.bai"
+ return "Results/STAR_Aln_2/" + wildcards.sample + "_rg_genomic_pp_rmdup_uniq_split.bam.bai"
# single end, no filter
elif len(sample_table["reverse_read"].dropna() ) == 0:
- return "Results/STAR_Aln_2/" + wildcards.sample + "_rg_genomic_rmdup_uniq_split.bai"
+ return "Results/STAR_Aln_2/" + wildcards.sample + "_rg_genomic_rmdup_uniq_split.bam.bai"
else :
raise Exception(wildcards.sample + " is sequenced in pair end & single end mode\n This workflow do not accept this case!\n")
@@ -39,8 +39,8 @@ rule ASEReadCounter:
input:
ref = config["fasta_ref"],
vcf = 'Results/hard_filtering/' + os.path.basename(config['vcf']).replace('.vcf.gz','') + '_SNP_filtered.vcf.gz',
- bam = get_splitNCigar_bam ,
- bai = get_splitNCigar_bai
+ bam = get_splitNCigar_bam,
+ bai = get_splitNCigar_bai
output:
cvs = 'Results/Counter/{sample}_ASEReadCounter.csv'
log:
@@ -52,6 +52,7 @@ rule ASEReadCounter:
baseQuality = config['baseQuality']
shell:
"""
+ ln -s {input.bam}.bai `sed 's/.bam.bai/.bai/'`
java -Xmx{params.mem} -jar {params.jar} -T ASEReadCounter -R {input.ref} -o {output} -I {input.bam} -sites {input.vcf} -U ALLOW_N_CIGAR_READS --minMappingQuality {params.mappingQuality} --minBaseQuality {params.baseQuality} 2> {log}
"""
@@ -61,7 +62,7 @@ rule phASER:
bam = get_splitNCigar_bam,
bai = get_splitNCigar_bai
output:
- cvs = 'Results/Counter/{sample}_phASER.csv'
+ csv = 'Results/Counter/{sample}_phASER.allelic_counts.txt'
log:
threads: config["phASER"]["cpu"]
params:
diff --git a/Snakemake/1000RNASeq_chicken/ASE/rules/basics.smk b/Snakemake/1000RNASeq_chicken/ASE/rules/basics.smk
index f37be08..18cc013 100644
--- a/Snakemake/1000RNASeq_chicken/ASE/rules/basics.smk
+++ b/Snakemake/1000RNASeq_chicken/ASE/rules/basics.smk
@@ -35,7 +35,7 @@ rule bam_index:
"{file}.bam"
output:
- temp("{file}.bai")
+ temp("{file}.bam.bai")
shell:
"samtools index {input} "
--
GitLab