SNAKE官网,snake snakes
用snakemake从零开始构建完整的甲基化信息分析流程。
看完文章的收获:专业的snakemake编程技巧;甲基化管道
作者拥有比利时根特大学的生物信息学硕士学位,拥有5年的分析师经验。他能迅速回应读者的问题,并提供技术支持。欢迎订阅专栏。
文章的前言,snakemake的介绍,甲基化生物信息的分析过程,snakemake代码的分发。
需求:在linux服务器上部署甲基化分析管道,用snakemake管理分析过程。这是基础部分,后面会带来详细的专业流程搭建方法。
一、snakemake简介snakemake起源于2012年。在其开发之前,市场上已经有很多开源的工作流引擎,比如Biopipe、Galaxy等。snakemake开发的目的是解决一些现有工作流引擎的缺点。最核心的好处是,只需要使用git就可以在远程服务器上部署盛鑫分析工作流,并以类似python语言语法的形式展现出来。有兴趣的话可以看看作者的文献。
Snakemake原创文学链接:点击阅读
二。甲基化生物信息分析流程甲基化信息分析的核心是基于illumina等测序技术仪器生成fastq数据文件,主要读取甲基化位点,即fastq文件中的碱基C,测序是核心。例子如下:
@
甲基化管道中文件处理频繁,用snakemake管理工作流方便易维护。让我们开始分发和编写代码并讨论它。
第三,分布式编写snakemake代码。接下来,一步一步写代码。
1.流程:fastp质量控制fastq文件:
Fastp可以输入一个fastq文件(单端测序)或两个fastq文件(双端测序)。
规则fastp:输入: c10ca-18r11313 _ r1.fastq.gz ,#输入第一个fastq文件 C10Ca-18R11313_R2.fastq.gz #输入第二个fastq文件输出: c10ca-18r11313 _ r1pre。#输出第一个fastq文件 c10ca-18r 11313 _ r2pre . fastq . gz #输出第二个fastq文件shell: mkdir-p test fastp-I { input[0]}-I { input[1]}-o { output[0]}-。
fastp-1c 10 ca-18r 11313 _ R1 . fastq . gz-1c 10 ca-18r 11313 _ R2 . fastq . gz-oc 10 ca-18r 11313 _ R1 pre . fastq . gz-oc 10 ca-18r 11313 _ R2 pre
构建工作的十克.使用shell:/usr/hxsdbbt/bash提供的核心:1要求更多线程的规则将被缩减。作业计数:计数作业1 fastp 1[周四10月14日17:23:16 2021]规则fastp:输入:C10Ca-18R11313_R1.fastq.gz,C10Ca-18R11313_R2.fastq.gz输出:C10Ca-18R11313_R1pre.fastq.gz,C10Ca-18r 11313 _ r2pre。fastq。gz jobid:0读取1过滤前:总读取数:98502总碱基数:14775300Q20碱基数:14338154(97.0414%)Q30碱基数:13494654(91.3325%)过滤后读取1:总读取数:97889总碱基数:14081251Q20碱基数:13702479(13702479-I C10Ca-18r 11313 _ R2。fastq。广州-O C10Ca-18r 11313 _ R1预。fastq。广州-O C10Ca-18r 11313 _ R2预。fastq。gz fastp v 0。14 .1,用时:3秒[2021年10月14日星期四17时23分19秒]已完成作业一步中的0.1步(100%)完成完整日志:/public/DUAN/methyl ization/.蛇的制造/日志/2021-10-14t 172416.90242435465蛇牌。原木流程2:俾斯麦软件比对,使用上面质控后的表达谱数据文件:
值得注意的是输出这一个规则里面没有,原因是俾斯麦这一步生成的是一个结果文件夹,包含结果文件,无法指定生成文件的名字,所以在这个规则里就省略了
规则俾斯麦:输入: C10Ca-18R11313_R1pre.fastq.gz , C10Ca-18r 11313 _ R2 pre。fastq。gz shell: /public/DUAN/甲基化/Bismark _ ref/Bismark-bow tie 2-path _ to _ bow tie/public/DUAN/甲基化/Bismark _ ref/bow tie 2-2。3 .4 .2-Linux-x86 _ 64-n0-l20-quiet-un-ambiguity-Sam-o/public/DUAN/methylation运行结果如下:
C在合格程序生成器(Certified Program Generator的缩写)上下文中甲基化:60.9%C在酒店集团上下文中甲基化:0.5%C在海南航空公司上下文中甲基化:0.5%C在未知上下文(CN或中国)中甲基化:0.1%俾斯麦在0d 0h 2m 15s内完成===============================俾斯麦运行完成=============================[2021年10月19日星期二11:09:16]已完成作业一步中的0.1步(100%)完成完成日志:/public/log snake make/log/2021-10-19t 110700.899050。蛇牌。原木流程3:重复数据删除_俾斯麦去除重复序列,使用上面比对生成的山姆文件:
规则Bismark _ dedup:输入: C10Ca-18r 11313 _ R1 pre _ Bismark _ bt2 _ PE。 Sam shell: /public/DUAN/metalization/duplicate _ Bismark-p { input }-output _ dir/public/DUAN/metalization 运行结果如下:
跳过标题行:@ SQ SN:chrX LN:155270560跳过标题行:@ SQ SN:chrY LN:59373566跳过标题行:@ SQ SN:chrMT LN:16569跳过标题行:@ PG ID:Bismark VN:v 0。20 .0 CL:俾斯麦-蝴蝶结2-path _ to _蝴蝶结/public/DUAN/甲基化/Bismark _ ref/领结2-2。3 .4 .2-Linux-x86 _ 64-Nsnake make/log/2021-10-19t 111433.699220。蛇牌。原木流程4:使用我修改过的封装好的大蟒文件,从山姆文件提取变异系数文件:
规则extract _ CpG:input: C10Ca-18r 11313 _ R1 pre _ bismark _ bt2 _ PE。重复。Sam shell: python/public/DUAN/DMR/extract _ CpG _ data。py-I { input }-o/public/DUAN/甲基化/C10Ca-18r 11313。封面 运行结果如下:
/public/DUAN/DMR/extract _ CpG _ data。py:51:弃用警告:“U”模式已弃用f=open(文件名, rU ){ }[2021年10月19日星期二11:22:00]已完成一个步骤中的0.1个(100%)已完成生成的变异系数文件格式如下:
后续就可以用决策支持系统(决策支持系统)读入变异系数文件进行分析组间甲基化差异啦