1000RNASeq ASE : add some changes from initial version uncommit

Snakemake/1000RNASeq_chicken/ASE/Snakefile

@@ -150,8 +150,12 @@ final_outputs.append('Results/Summary/' + os.path.basename(config['vcf']).replac
 
 # counter
 for sample in table["sample_name"].unique():
-	final_outputs.append('Results/Counter/' + sample + '_ASEReadCounter.csv')
-	final_outputs.append('Results/Counter/' + sample + '_phASER.allelic_counts.txt')
+	final_outputs.append('Results/ASEReadCounter/' + sample + '_ASEReadCounter.csv')
+	final_outputs.append('Results/phASER/' + sample + '_phASER.allelic_counts.txt')
+	final_outputs.append('Results/phASER_gene_ae/' + sample + '_phASER.gene_ae.txt')
+
+final_outputs.append('Results/phASER_pop/' + os.path.basename(config['vcf']).replace('.vcf.gz','_phASER.bed'))
+final_outputs.append('Results/phASER_pop/' + os.path.basename(config['vcf']).replace('.vcf.gz','_phASER.gw_phased.bed'))
 
 # print(final_outputs)
 # ['Results/STAR_Aln_2/samplePE_rg_genomic_pp_rmdup_uniq_split.bam',

Snakemake/1000RNASeq_chicken/ASE/resources_SLURM.yaml

@@ -37,5 +37,8 @@ SplitNCigarReads :
 phASER:
     mem : "20G"
 
+phASER_expr_matrix:
+    cpu : 8
+
 ASEReadCounter:
     mem : "20G"

Snakemake/1000RNASeq_chicken/ASE/rules/ASE_counter.smk

@@ -44,7 +44,7 @@ rule ASEReadCounter:
 		bam = get_splitNCigar_bam,
 		bai = get_splitNCigar_bai
 	output:
-		cvs = 'Results/Counter/{sample}_ASEReadCounter.csv'
+		cvs = 'Results/ASEReadCounter/{sample}_ASEReadCounter.csv'
 	params:
 		mem  = config["ASEReadCounter"]["mem"],
 		mappingQuality = config['mappingQuality'],
@@ -54,6 +54,17 @@ rule ASEReadCounter:
 		gatk --java-options "-Xmx{params.mem}" ASEReadCounter -R {input.ref} -O {output} -I {input.bam} --variant {input.vcf}   --min-mapping-quality {params.mappingQuality} --min-base-quality {params.baseQuality}
 		"""
 
+rule gtf2bed:
+	input:
+		gtf = config['gtf_ref']
+	output:
+		bed = 'Results/ref_bed/' + os.path.basename(config['gtf_ref']).replace('.gtf', '.bed')
+	script:
+		"""
+		../script/gtf2bed.py
+		"""
+
+
 rule phASER:
 	input:
 		vcf = 'Results/vcfFiltration/' + os.path.basename(config['vcf']).replace('.vcf.gz','_FsQdBiall.vcf.gz'),
@@ -61,7 +72,8 @@ rule phASER:
 		bam = get_splitNCigar_bam,
 		bai = get_splitNCigar_bai
 	output:
-		csv = 'Results/Counter/{sample}_phASER.allelic_counts.txt'
+		csv = 'Results/phASER/{sample}_phASER.allelic_counts.txt',
+		hap_count = 'Results/phASER/{sample}_phASER.haplotypic_counts.txt'
 	threads:  config["phASER"]["cpu"]
 	params:
 		mappingQuality = config['mappingQuality'],
@@ -72,3 +84,32 @@ rule phASER:
 		"""
 		python2.7 {config[bin_dir]}/phaser.py --vcf {input.vcf} --bam {input.bam} --sample {wildcards.sample} --threads {threads} --o Results/Counter/{wildcards.sample}_phASER --paired_end {params.paired_option} --baseq {params.baseQuality} --mapq {params.mappingQuality} --id_separator {params.idSeparator}
 		"""
+
+rule phASER_gene_ae:
+	input:
+		hap_count = "Results/phASER/{sample}_phASER.haplotypic_counts.txt",
+		bed = 'Results/ref_bed/' + os.path.basename(config['gtf_ref']).replace('.gtf', '.bed')
+	output:
+		gene_ae = "Results/phASER_gene_ae/{sample}_phASER.gene_ae.txt"
+	params:
+		idSeparator = '_' if config['id_separator'] == '' else config['id_separator']
+	shell:
+		"""
+		python2.7 {config[bin_dir]}/phaser_gene_ae.py --haplotypic_counts {input.hap_count} --features {input.bed} --o {output} --id_separa
+		"""
+
+rule phASER_expr_matrix:
+	input:
+		expand("Results/phASER_gene_ae/{sample}_phASER.gene_ae.txt", sample=table["sample_name"].unique()),
+		bed = 'Results/ref_bed/' + os.path.basename(config['gtf_ref']).replace('.gtf', '.bed')
+	output:
+		'Results/phASER_pop/' + os.path.basename(config['vcf']).replace('.vcf.gz','_phASER.bed'),
+		'Results/phASER_pop/' + os.path.basename(config['vcf']).replace('.vcf.gz','_phASER.gw_phased.bed')
+	params:
+		gene_ae_dir = "Results/phASER_gene_ae/",
+		out_dir = 'Results/phASER_pop/'
+	threads: config["phASER_expr_matrix"]["cpu"]
+	shell:
+		"""
+		python2.7 {config[bin_dir]}/phaser_expr_matrix.py --t {threads} --gene_ae_dir {params.gene_ae_dir} --feature {input.bed} --o {param
+		"""