上一次讲述了如何利用snakemake进行RNA-seq分析流程的构建：
RNAseq分析流程：从测序Rawdata到Count输出，但是只针对单端测序数据。这次我又改进了一点，对双端测序结果也适用。
cat pairEndsnakemake
configfile: "pairedEnd.yaml"
rule all:
input:
expand("results/{sample}.count",sample=config["samples"])
rule bwa:
input:
fw=lambda wildcards: expand("samples/{sample}_R1.fq.gz", sample=wildcards.sample),
re=lambda wildcards: expand("samples/{sample}_R2.fq.gz", sample=wildcards.sample)
params:INDEX=config["annotation"]["reference"]
output:
"results/{sample}.sam"
shell:
"bowtie2 -x {params.INDEX} -1 {input.fw} -2 {input.re} -S {output}"
rule sam2bam:
input:
"results/{sample}.sam"
output:
"results/{sample}.bam"
shell:
"samtools view -Sb {input} > {output}"
rule samSort:
input:
"results/{sample}.bam"
output:
"sorted/{sample}_sort.bam"
shell:
"samtools sort {input} -o {output}"
rule bamIndex:
input:
"sorted/{sample}_sort.bam"
output:
"results/{sample}.bam.bai"
shell:
"samtools index {input} {output}"
rule bamCount:
input:
sample="sorted/{sample}_sort.bam",
annotation=config["annotation"]["gtf"]
output:
"results/{sample}.count"
shell:
"htseq-count -f bam -r name -s no {input.sample} {input.annotation} > {output}"
有需要代码指导的可以联系我，我也会持续更新，先暂时讲到这里。
