Skip to content
Snippets Groups Projects
Commit c8a04e87 authored by Shaman Narayanasamy's avatar Shaman Narayanasamy
Browse files

Add Analysis MG.rules

parent eb358b43
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 406 additions and 0 deletions
rules/Analysis/MG.rules 0 → 100644
+ 406
− 0
View file @ c8a04e87
rule ANALYSIS_ANNOTATE:
log:
AN_LOG
benchmark:
"%s/benchmarks/ANALYSIS_ANNOTATE.json" % AN_OUT
input:
'{dir}/MG.assembly.merged.fa'.format(dir=A_OUT),
expand("{path}/{db}", path=DBPATH, db=config["prokka"]["databases"])
output:
"%s/annotation/annotation.filt.gff" % AN_OUT
shell:
"""
prokka --force --outdir {AN_OUT}/annotation --cpus {THREADS} --metagenome --norrna {input[0]} >> {log} 2>&1
# Prokka gives a weird gff file with all the sequences. We need to write some small code to produce a file that
# cleans up the output
LN=`grep -Hn "^>" {AN_OUT}/annotation/PROKKA_*.gff | head -n1 | cut -f2 -d ":"`
LN1=1
LN=$(($LN-$LN1))
head -n $LN {AN_OUT}/annotation/*.gff | grep -v "^#" | sort | uniq | grep -v "^==" > {output}
"""
rule ANALYSIS_MG_GENE_COVERAGE:
log:
AN_LOG
benchmark:
"%s/benchmarks/ANALYSIS_MG_GENE_COVERAGE.json" % AN_OUT
input:
"%s/annotation/annotation.filt.gff" % AN_OUT,
"%s/MG.reads.sorted.bam" % A_OUT
output:
"%s/MG.annotation.bed" % AN_OUT,
"%s/MG.gene_depth.hist" % AN_OUT,
"%s/MG.gene_depth.avg" % AN_OUT,
"%s/MG.gene.len" % AN_OUT,
"%s/MG.prokkaID2ec.txt" % AN_OUT
shell:
"""
coverageBed -hist -abam {input[1]} -b {input[0]} | grep -v "^all" > {output[0]}
paste <(cat {output[0]} | cut -f9 | cut -f1 -d \";\" | sed -e \"s/ID=//g\") \
<(cut -f10,11,12,13 {output[0]}) > {output[1]}
## This code was adapted and modified from the CONCOCT script to calculate depth
## It prints out a file that contains the average depth of all the genes
awk -v OFS='\t' 'BEGIN {{pc=""}}
{{
c=$1;
if (c == pc) {{
cov=cov+$2*$5;
}} else {{
print pc,cov;
cov=$2*$5;
pc=c
}}
}} END {{print pc,cov}}' < {output[1]} | tail -n +2 > {output[2]}
# Record gene length file
cut -f 1,4 {output[1]} | uniq > {output[3]}
grep "eC_number=" {input[0]} | cut -f9 | cut -f1,2 -d ';'| sed 's/ID=//g'| sed 's/;eC_number=/\t/g' > {output[4]}
"""
rule ANALYSIS_MG_VARIANT_CALLING:
log:
AN_LOG
benchmark:
"%s/benchmarks/ANALYSIS_MG_VARIANT_CALLING.json" % AN_OUT
input:
"%s/MG.assembly.merged.fa" % A_OUT,
"%s/MG.reads.sorted.bam" % A_OUT,
output:
"%s/MG.variants.isec.vcf.gz" % AN_OUT
shell:
"""
echo "[x] VARIANT CALLING `date +"%Y/%m/%d %H:%M:%S"`"
if [[ ! -f {input[1]}.bai ]]
then
echo "Bam index doesn't exist, Creating one..."
echo "Indexing bam: {input[1]}"
samtools index {input[1]}
fi
if [[ ! -f {input[0]}.fai ]]
then
echo "Fasta index doesn't exist, Creating one..."
echo "Indexing fasta: {input[0]}"
samtools faidx {input[0]}
fi
#temporary directory and files
VCF_MPU=$(mktemp --tmpdir={TMPDIR} -t "XXXXXX.mpu.vcf")
VCF_FRB=$(mktemp --tmpdir={TMPDIR} -t "XXXXXX.frb.vcf")
VCF_PLT=$(mktemp --tmpdir={TMPDIR} -t "XXXXXX.plt.vcf")
### run_mpileup {input[0]} {input[1]} $VCF_MPU
samtools mpileup -uf {input[0]} {input[1]} | bcftools view -vcg - > $VCF_MPU
bgzip -c $VCF_MPU > $VCF_MPU.gz
tabix -f -p vcf $VCF_MPU.gz
### run_freebayes {input[0]} {input[1]} $VCF_FRB
freebayes -f {input[0]} {input[1]} > $VCF_FRB
bgzip -c $VCF_FRB > $VCF_FRB.gz
tabix -f -p vcf $VCF_FRB.gz
### run_platypus {input[0]} {input[1]} $VCF_PLT
Platypus.py callVariants --refFile={input[0]} --bamFiles={input[1]} --nCPU={THREADS} -o $VCF_PLT
bgzip -c $VCF_PLT > $VCF_PLT.gz
tabix -f -p vcf $VCF_PLT.gz
### "Merging outputs "
## Must remove colons from the contig names in upstream steps. Unable to merge the variants
## due to this problem
vcf-isec -f -a -n +2 $VCF_PLT.gz $VCF_FRB.gz $VCF_MPU.gz > {AN_OUT}/MG.variants.isec.vcf
# Compress and index the output.
bgzip -c {AN_OUT}/MG.variants.isec.vcf > {AN_OUT}/MG.variants.isec.vcf.gz
tabix -f -p vcf {AN_OUT}/MG.variants.isec.vcf.gz
# Clean up directory
echo "Cleaning up directory"
cat log.txt >> {log}
rm -f {AN_OUT}/MG.variants.isec.vcf log.txt
"""
rule ANALYSIS_MG_CONTIG_LENGTH:
log:
AN_LOG
benchmark:
"%s/benchmarks/ANALYSIS_MG_CONTIG_LENGTH.json" % AN_OUT
input:
"%s/MG.assembly.merged.fa" % A_OUT,
output:
"%s/MG.assembly.length.txt" % AN_OUT,
"%s/MG.assembly.gc_content.txt" % AN_OUT,
shell:
"""
echo "[x] LENGTH `date +"%Y/%m/%d %H:%M:%S"`" >> {log}
echo "Obtaining contig lengths"
perl {SRCDIR}/fastaNamesSizes.pl {input} > {output[0]}
echo "Obtaining GC content"
TMP_GC=$(mktemp --tmpdir={TMPDIR} -t "gc_out_XXXXXX.txt")
perl {SRCDIR}/get_GC_content.pl {input} $TMP_GC
# Th program above provides a file gc_out.txt. This command cleans the output
echo "Clean up output"
cut -f1,2 $TMP_GC | sed -e 's/>//g'> {output[1]}
echo "Remove intermediate files"
rm $TMP_GC
"""
rule ANALYSIS_MG_CONTIG_COVERAGE:
log:
AN_LOG
benchmark:
"%s/benchmarks/ANALYSIS_MG_CONTIG_COVERAGE.json" % AN_OUT
input:
"%s/MG.reads.sorted.bam" % A_OUT,
"%s/MG.assembly.merged.fa" % A_OUT,
output:
"%s/MG.assembly.contig_coverage.txt" % AN_OUT,
"%s/MG.assembly.contig_depth.txt" % AN_OUT,
"%s/MG.assembly.contig_flagstat.txt" % AN_OUT,
shell:
"""
echo "[x] COVERAGE AND DEPTH `date +"%Y/%m/%d %H:%M:%S"`" >> {log}
echo "Creating genome file ..." >> {log}
if [[ ! -f {input[1]}.fai ]]
then
echo "No fasta index! Creating one." >> {log}
samtools faidx {input[1]}
fi
cat {input[1]}.fai | awk '{{print $1 \"\t0\t\" $2}}' > {input[1]}.bed3
echo "Done creating bed file" >> {log}
echo "Running BEDTools coverage calculation ..." >> {log}
coverageBed -abam {input[0]} -b {input[1]}.bed3 > {output[0]}
echo "Coverage calculation done" >> {log}
echo "Running BEDTools for average depth in each position" >> {log}
TMP_DEPTH=$(mktemp --tmpdir={TMPDIR} -t "depth_file_XXXXXX.txt")
genomeCoverageBed -ibam {input[0]} | grep -v "genome" > $TMP_DEPTH
echo "Depth calculation done" >> {log}
## This method of depth calculation was adapted and modified from the CONCOCT code
awk -v OFS='\t' 'BEGIN {{pc=""}}
{{
c=$1;
if (c == pc) {{
cov=cov+$2*$5;
}} else {{
print pc,cov;
cov=$2*$5;
pc=c}}
}} END {{print pc,cov}}' $TMP_DEPTH | tail -n +2 > {output[1]}
echo "Remove the temporary file" >> {log}
rm $TMP_DEPTH
echo "flagstat" >> {log}
samtools flagstat {input[0]} | cut -f1 -d ' ' > {output[2]}
"""
def mg_fqfiles():
raw = expand('{dir}/{raw}', raw=['MG.R1.fq', 'MG.R2.fq'], dir=P_OUT)
uniq = expand('{dir}/{uniq}', uniq=['MG.R1.uniq.fq', 'MG.R2.uniq.fq'], dir=P_OUT)
trim = expand('{dir}/{trim}', trim=[
'MG.R1.uniq.trimmed.fq',
'MG.R2.uniq.trimmed.fq',
'MG.SE.uniq.trimmed.fq'], dir=P_OUT)
filter = expand('{dir}/{filter}', filter=expand([
'MG.R1.uniq.trimmed.{f}.fq',
'MG.R2.uniq.trimmed.{f}.fq',
'MG.SE.uniq.trimmed.{f}.fq'], f=config['human_filtering']['filter']), dir=P_OUT)
if config['preprocessing_filtering']:
return raw + uniq + trim + filter
return raw + uniq + trim
rule ANALYSIS_MG_READ_COUNT:
log:
AN_LOG
benchmark:
"%s/benchmarks/ANALYSIS_MG_READ_COUNT.json" % AN_OUT
input:
mg_fqfiles()
output:
'%s/MG.read_counts.txt' % AN_OUT
run:
for idx, f in enumerate({input}):
if idx == 0:
shell("wc -l {f} > {output}")
else:
shell("wc -l {f} >> {output}")
rule ANALYSIS_VIZBIN:
log:
AN_LOG
benchmark:
"%s/benchmarks/ANALYSIS_VIZBIN.json" % AN_OUT
input:
"%s/MG.assembly.merged.fa" % A_OUT,
output:
"%s/MG.vizbin.points" % AN_OUT,
"%s/MG.vizbin.filtered.fa" % AN_OUT,
"%s/MG.vizbin.with-contig-names.points" % AN_OUT
shell:
"""
TMP_VIZBIN=$(mktemp --tmpdir={TMPDIR} -dt "VIZBIN_XXXXXX")
# tools used in this script: they must be in the PATH
VIZBIN={SRCDIR}/vizbin.py
echo "[x] VIZBIN `date +"%Y/%m/%d %H:%M:%S"`" >> {log}
python $VIZBIN {input} {AN_OUT}/MG.vizbin {config[vizbin][dimension]} \
-l {config[vizbin][length]} \
-s {config[vizbin][size]} \
-t {config[vizbin][theta]} \
-p {config[vizbin][perp]} \
-f {config[vizbin][minlen]} \
-o {output[1]} >> {log}
cp {AN_OUT}/MG.vizbin.{config[vizbin][dimension]}.dat $TMP_VIZBIN/data.dat
BACK=$PWD
cd $TMP_VIZBIN
bh_tsne >> $BACK/{log}
cd $BACK
cp $TMP_VIZBIN/points.txt {AN_OUT}/MG.vizbin.points
rm -rf $TMP_VIZBIN
grep '^>' {output[1]} | sed 's/>//g' > {AN_OUT}/MG.vizbin.contigs
paste {AN_OUT}/MG.vizbin.contigs {output[0]} | sed -e 's/,/\t/g' > {output[2]}
"""
def analysis_plot_files_output():
return expand('{dir}/{name}', name=[
"IMP-reads_density.png",
"IMP-rpkm_density.png",
"IMP-coverage_density.png",
"IMP-depth_density.png",
"IMP-var_count.png",
"IMP-var_density.png",
"IMP-vizbin_length.png",
"IMP-vizbin_length_GC.png",
"IMP-vizbin_length_MGcov.png",
"IMP-vizbin_length_MGdepth.png",
"IMP-vizbin_length_MGvardens.png",
"IMP-vizbin_length_geneDensity.png",
"IMP-vizbin_standard.png",
"MG.read_stats.html",
"MG.read_stats.txt",
"MG_mapping_stats.html",
"MG_mapping_stats.txt",
"assembly_stats.html",
"assembly_stats.txt"], dir='%s/results' % AN_OUT)
def analysis_stats_files_output():
return expand([
'{dir}/{stat_flag}/{rtype}/cycle_composition_{n}.{ext}',
'{dir}/{stat_flag}/{rtype}/cycle_quality_{n}.{ext}',
'{dir}/{stat_flag}/{rtype}/cycle_quality_box_{n}.{ext}',
'{dir}/{stat_flag}/{rtype}/info.tab',
'{dir}/{stat_flag}/{rtype}/lane_tile_quality_{n}.{ext}',
'{dir}/{stat_flag}/{rtype}/quality_QQ.{ext}',
'{dir}/{stat_flag}/{rtype}/reads_length.{ext}',
'{dir}/{stat_flag}/{rtype}/reads_quality.{ext}'],
n = ['1', '2'],
ext = ['gnuplot', 'png', 'tab'],
dir = P_OUT,
stat_flag = ['stats', 'stats_after_preprocessing'],
rtype = ['MG'])
rule ANALYSIS_PLOT:
log:
AN_LOG
benchmark:
"%s/benchmarks/ANALYSIS_PLOT.json" % AN_OUT
input:
expand('{dir}/{name}', name=[
'MG.read_counts.txt',
'MG.assembly.contig_flagstat.txt',
'MG.assembly.contig_coverage.txt',
'MG.assembly.contig_depth.txt',
'MG.variants.isec.vcf.gz',
'MG.assembly.gc_content.txt',
"MG.vizbin.with-contig-names.points",
"annotation/annotation.filt.gff"], dir=AN_OUT)
output:
analysis_plot_files_output()
params:
outdir = "%s/results" % AN_OUT
shell:
"""
PLOT_SCRIPT="{SRCDIR}/IMP_visualize_MG.R"
echo "[x] PLOT `date +"%Y/%m/%d %H:%M:%S"`" >> {log}
mkdir -p {AN_OUT}/results
Rscript $PLOT_SCRIPT {AN_OUT}/results {input}
"""
rule ANALYSIS_KRONA_PLOT_MG:
log:
AN_LOG
benchmark:
"%s/benchmarks/ANALYSIS_KRONA_PLOT.json" % AN_OUT
input:
"%s/MG.prokkaID2ec.txt" % AN_OUT,
"%s/ec2pwy.txt" % U_OUT,
"%s/pwy2hierarchy.txt" % U_OUT,
"%s/MG.gene_depth.avg" % AN_OUT,
"%s/MG.gene.len" % AN_OUT
output:
"%s/results/MG.gene_kegg_krona.txt" % AN_OUT,
"%s/results/MG.gene_kegg_krona.html" % AN_OUT
shell:
"""
echo "[x] PLOT KRONA `date +"%Y/%m/%d %H:%M:%S"`" >> {log}
echo {input}
echo {output[0]}
python {SRCDIR}/genes.to.kronaTable.py -i {input[0]} -H {input[1]} -m {input[2]} -c {input[3]} -L {input[4]} -o {output[0]}
ktImportText -o {output[1]} {output[0]}
"""
rule ANALYSIS_MG_QUALITY_STATS:
input:
expand('{dir}/{raw}', raw=['MG.R1.fq', 'MG.R2.fq'], dir=P_OUT)
output:
expand(['{dir}/cycle_composition_{n}.{ext}',
'{dir}/cycle_quality_{n}.{ext}',
'{dir}/cycle_quality_box_{n}.{ext}',
'{dir}/info.tab',
'{dir}/lane_tile_quality_{n}.{ext}',
'{dir}/quality_QQ.{ext}',
'{dir}/reads_length.{ext}',
'{dir}/reads_quality.{ext}'
], n=['1', '2'], ext=['gnuplot', 'png', 'tab'], dir='%s/stats/MG' % P_OUT)
benchmark:
"%s/benchmarks/PREPROCESSING_MG_QUALITY_STATS.json" % P_OUT
log:
P_LOG
shell:
"""
ht2-stat --encode=sanger -q -P -t {THREADS} -o {P_OUT}/stats/MG -i {input} >> {log} 2>&1
ht2-stat-draw.pl --dir {P_OUT}/stats/MG >> {log} 2>&1
"""
rule ANALYSIS_MG_PREPROCESSED_QUALITY_STATS:
log:
P_LOG
input:
expand('{dir}/{trim}', trim=[
'MG.R1.uniq.trimmed.fq',
'MG.R2.uniq.trimmed.fq'], dir=P_OUT)
output:
expand(['{dir}/cycle_composition_{n}.{ext}',
'{dir}/cycle_quality_{n}.{ext}',
'{dir}/cycle_quality_box_{n}.{ext}',
'{dir}/info.tab', '{dir}/lane_tile_quality_{n}.{ext}',
'{dir}/quality_QQ.{ext}',
'{dir}/reads_length.{ext}',
'{dir}/reads_quality.{ext}'
], n=['1', '2'], ext=['gnuplot', 'png', 'tab'], dir='%s/stats_after_preprocessing/MG' % P_OUT)
shell:
"""
ht2-stat --encode=sanger -q -P -t {THREADS} -o {P_OUT}/stats_after_preprocessing/MG -i {input} >> {log} 2>&1
ht2-stat-draw.pl --dir {P_OUT}/stats_after_preprocessing/MG >> {log} 2>&1
"""
benchmark:
"%s/benchmarks/PREPROCESSING_MG_PREPROCESSED_QUALITY_STATS.json" % P_OUT
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