一次snakemake分析数据的经验

参考文章：
snakemake官方教程文档
https://blog.csdn.net/u012110870/article/details/85330457
https://www.jianshu.com/p/14b9eccc0c0e
www.php-master.com/post/352109.html
其中第一个和第三个基本一样，百度一下基本上就是这些例子，按照这些例子运行一下基本上都没问题，我就又从这里下载了RNA-seq数据（用snakemake写RNA-Seq流程
），然后自己在上面的基础上，编写了一个简单的snakemake流程，说是简单，也折腾了2-3天，主要是配置文件搞不对。
这是下载的数据：

ls samples/*.gz
samples/control_R1.fq.gz  samples/treated_R1.fq.gz
samples/control_R2.fq.gz  samples/treated_R2.fq.gz

我打算进行这样的分析流程，其实很简单：比对--转换bam文件--排序--建索引

分析流程

这是输入的配置文件，主要是文件的名称和路径

cat config.yaml #关于配置文件的介绍，网上一大把资源，重点就是格式，左对齐，以空格分隔
samples:
 control_R1: samples/control_R1.fq.gz
 control_R2: samples/control_R2.fq.gz
 treated_R1: samples/treated_R1.fq.gz
 treated_R2: samples/treated_R2.fq.gz

其次就是分析的开始了，这是我当前的目录文件：

 tree reference/ samples/
reference/#根据参考文章建立的19号染色体索引
├── chr19.1.bt2
├── chr19.2.bt2
├── chr19.3.bt2
├── chr19.4.bt2
├── chr19.fa
├── chr19.fa_bowtie2-build.log
├── chr19.rev.1.bt2
├── chr19.rev.2.bt2
├── GRCm38.83.chr19.gtf
└── nohup.out
samples/#下载的测序数据
├── control_R1.fq.gz
├── control_R2.fq.gz
├── nohup.out
├── treated_R1.fq.gz
├── treated_R2.fq.gz
├── wget-log
└── wget-log.1

下面就开始编写分析流程

rule all:
        input:#这里只要通过expand告诉分析流程你最终输出文件{sample}对应的什么就可以了
              #所以我这里的前3行都是多余的，我是一行行试，才发现多余的 ○|~|_
        #       expand("results/{sample}.sam",sample=config["samples"]), 
        #       expand("results/{sample}.bam",sample=config["samples"]), 
        #       expand("sorted/{sample}_sort.bam",sample=config["samples"])
                  expand("results/{sample}.bam.bai",sample=config["samples"])
rule bwa:
        input:
                lambda wildcards: config["samples"][wildcards.sample]#伪函数，利用wildcards.sample由注释文件中获得sample对应的路径

        params:
                INDEX="reference/chr19"#参考19号染色体索引路径，以参数形式录入分析流程

        output:
                "results/{sample}.sam" #定义输出
        shell:
                "bowtie2 -x {params.INDEX} -q {input} -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}"

最终分析的结果在sorted，和results两个文件中，还可以通过temp，protect函数进行删除和保护，看一下两个文件夹里有些什么：

 tree results/ sorted/
results/
├── control_R1.bam
├── control_R1.bam.bai
├── control_R1.sam
├── control_R2.bam
├── control_R2.bam.bai
├── control_R2.sam
├── treated_R1.bam
├── treated_R1.bam.bai
├── treated_R1.sam
├── treated_R2.bam
├── treated_R2.bam.bai
└── treated_R2.sam
sorted/
├── control_R1_sort.bam
├── control_R2_sort.bam
├── treated_R1_sort.bam
└── treated_R2_sort.bam

在写流程的时候，还有个小trick，可以通过snakemake -s Snakefile -np，跑一下各个rule，看一下正确，再除去-np，重新开始跑程序，会很省时省力。

但是，使用shell脚本也可以达到类似的目的

 ls *.gz|while read id; do bowtie2 -x ../reference/chr19 -q ${id} |samtools view -Sb -| samtools sort |samtools index - ${id}.bai ;done

 ls -lh samples/*.bai #不过只能得到最终的输出文件.bai，中间的过程文件不知道怎么获得
-rw-r--r-- 1 root root 55K May 15 09:28 samples/control_R1.fq.gz.bai
-rw-r--r-- 1 root root 55K May 15 09:29 samples/control_R2.fq.gz.bai
-rw-r--r-- 1 root root 56K May 15 09:30 samples/treated_R1.fq.gz.bai
-rw-r--r-- 1 root root 56K May 15 09:31 samples/treated_R2.fq.gz.bai