RNAseq分析流程（snakemake+ script)

接之前的文章
RNAseq分析流程：从测序Rawdata到Count输出
一次snakemake分析数据的经验
官方文档
用snakemake写分析流程，其中一个优点可能就是可以插入pyhon.py和rscript.R
external-scripts

1.插入python 的script

将snakefile中的merge写成python脚本

import os
  import pandas as pd
  dm=pd.DataFrame()
  for filename in snakemake1.input:#遍历文件夹下的数据文件
    df=pd.read_csv(filename,sep="\t",names=["gene",os.path.basename(filename)])
    if(any(dm)):
        dm=pd.merge(dm,df,on="gene")
    else:
        dm=df
    
  dm.to_csv(snakemake1.output[0],sep="\t")#输出到文件夹下，并在ouput中定义好参数

在snakefile中 rule all替换成这样

rule all:
 input:
 #告诉snake我处理的是文件夹下的一系列文件
  expand("count/{sample}.count",sample=config["samples"]) output:
  "count/merge.count"
 script:
  "merge.py"

然后snakemake运行即可

2.插入R 的script

将snakefile中的merge写成R脚本

dm=c()
for(filename in snakemake@input){
df=read.csv(filename,sep="\t",col.names=c("gene",basename(filename)))
if(is.null(dm)){dm=df}
else{dm=merge(dm,df,by="gene")}
}

require(Desq2)
dm=dm[-c(1:5),]
write.table(dm,snakemake@output[[1]],sep="\t",quote=F,row.names=F)
#这里snakemake@output[[1]]不同于python的使用方式

在snakefile中将rule all替换成这样

rule all:
 input:
  expand("count/{sample}.count",sample=config["samples"])
 output:
  "count/merge.count"
 script:
  "merge.R"

均可以在count/文件夹下产生merge.count文件