参考文章：
snakemake官方教程文档
https://blog.csdn.net/u012110870/article/details/85330457
//www.greatytc.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
