rule strand_specific:
input:
bed = hg38_bed12,
bam = op.join(base, gse,"results/mapping", "{sample}", "{sample}.rmdup.bam"),
output:
op.join(base, gse,"results/strand_specific", "{sample}.strand_specific.txt")
params:
op.join(base, gse,"results/strand_specific")
shell:"""
mkdir -p {params}
echo "sample name:{wildcards.sample}" > {output}
infer_experiment.py -r {input.bed} -i {input.bam} >> {output} 2>&1
"""
rule bam2bw:
input:
bam = op.join(base, gse,"results/mapping", "{sample}", "{sample}.rmdup.bam"),
strand = op.join(base, gse,"results/strand_specific", "{sample}.strand_specific.txt")
output:
bw = op.join(base, gse, "results/bw", "{sample}.cpm.bw"),
fwd = op.join(base, gse, "results/bw", "{sample}.cpm.fwd.bw"),
rev = op.join(base, gse, "results/bw", "{sample}.cpm.rev.bw"),
threads:
5
log:
op.join(base, gse, "results/bw", "{sample}.log"),
shell:"""
strand1=$(grep -Po '(?<=":\s).+' {input.strand} | xargs | awk '{{if($1/($2+0.001)>2)print "reverse"; else print "forward"}}')
if [ strand1=="forward" ];then strand2="reverse";else strand2="forward";fi
bamCoverage -p {threads} -b {input.bam} --normalizeUsing CPM -bs 1 -o {output.bw} > {log} 2>&1
bamCoverage -p {threads} -b {input.bam} --normalizeUsing CPM -bs 1 --filterRNAstrand $strand1 -o {output.fwd} >> {log} 2>&1
bamCoverage -p {threads} -b {input.bam} --normalizeUsing CPM -bs 1 --filterRNAstrand $strand2 -o {output.rev} >> {log} 2>&1
"""
bam文件转bw文件并区分正负链
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。
推荐阅读更多精彩内容
- 整理ChIP-seq / CUT & Tag 分析时用到的工具。本文只对使用的工具用法进行简单介绍。 deepto...
- bedtools工具包bamtofastq转换工具 picard中的SamToFastq转换工具 samtools...
- 需求: PacBio ccs测序下机数据是bam格式的。而组装软件一般都需要fasta格式,因此需要将源文件转换成...
- 需求 在实际工作中,我经常要处理其他同事发来的csv文件。这些文件不是一个两个,而是十几二十个。因此,怎么准确区分...