bam文件转bw文件并区分正负链


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
    """
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容