搜索
查看: 12960|回复: 23

[WGBS] BS-seq 分析流程初学1

[复制链接]

1

主题

10

帖子

257

积分

中级会员

Rank: 3Rank: 3

积分
257
发表于 2017-2-24 00:27:09 | 显示全部楼层 |阅读模式
刚开始学习分析BS-seq数据,自己通过看文献,探索了一个粗浅的流程,将逐步与大家分享交流。希望能够互相交流与等到指导。
我分析的物种是斑马鱼(danio rerio)
1.所需软件:trimmomatic(预处理、reads质量控制)、bismark(BS-seq mapping)、bsseq (R package,Find DMRs)、ChipPeakanno(annotation R package)。
2.trimmomatic预处理,非常常规的参数:
[Shell] 纯文本查看 复制代码
for i in `ls *1.fastq.gz`
do x=${i/1.fastq.gz/}
echo java -jar /home/wenhao/work/biosoft/Trimmomatic-0.33/trimmomatic-0.33.jar PE -threads 12 -phred33 \
${x}"1.fastq.gz" ${x}"2.fastq.gz" ${x}"1.clean.fastq.gz" ${x}"1.unpaired.fastq.gz" ${x}"2.clean.fastq.gz" \
${x}"2.unpaired.fastq.gz" ILLUMINACLIP:/opt/biosoft/Trimmomatic-0.33/adapters/TruSeq3-PE.fa:2:30:10 \
AVGQUAL:20 TRAILING:20 MINLEN:50
done > trimmomatic.sh

3.bismark
3.1准备基因组文件:
[AppleScript] 纯文本查看 复制代码
bismark_genome_preparation --path_to_bowtie /opt/biosoft/bowtie2-2.2.9 --verbose --bowtie2 /home/data/genome_file_dir
#--path_to_bowtie 指定bowtie所在的路径
#--bowtie2使用bowtie2
#/home/data/genome_file_dir 参考基因组所在的路径


3.2
[AppleScript] 纯文本查看 复制代码
bismark --bowtie2 -p 30 --score_min L,0,-0.6 --phred33-quals -N 1 /home/wenhao/BS-seq-temperature-muscle -1 D2-TE22_ATCACG_L008_R1.clean.trim.fastq.gz -2 D2-TE22_ATCACG_L008_R2.clean.trim.fastq.gz
# -p 线程数
#--score_min ,bismark调用bowtie2,默认的参数是L,0,-0.6,但是这个时候map efficiency只有45%,低到离谱。查看bowtie2的默认参数是L,-0.6,-0.6,猜出bismark的比对过程可能过于严格,就将参数设置到了L,0,-0.6,这时候map efficiency升到的65%。虽然不是太高,但是已经提高不少。
#-N seed alignment允许的错配个数。一般设置为1

3.3 从bam文件中提取甲基化结果
[AppleScript] 纯文本查看 复制代码
for i in `ls *.bam`
do
echo bismark_methylation_extractor \
--paired-end \                     #成对测序
--no_overlap \                    #将read1和read2重合的部分剔除,不计算在内,防止出现偏差
--ignore 10 \                       # read1 5'端10bp,不计算在内。如果前10bp的ATGC比值分布比较怪异,需要将其除掉
--ignore_r2 10 \                  # read2 5'端10bp,不计算在内
--comprehensive \              
--merge_non_CpG \
--gzip \
--multicore 10 \
--bedGraph \          
--cytosine_report \
--scaffolds \                          #如果基因组有很多小的contig或者scaffold,这个时候必须加上,不然容易报错。斑马鱼这种模式生物还需要加上,我想大部分的基因组都需要加上,但是如果加上这个产生,将会是时间大大延长。最好的方法是bismark mapping比对之前就将基因组中的小片段给去除掉,不考虑在内
--counts \
--buffer_size 60G \
--genome_folder /home/data/genome_file_dir \
$i
done > bismark_methylation_extractor.sh


##未完待续 2017/2/24

























上一篇:2008年就有68种做富集分析的工具了!
下一篇:生信编程直播第8题-几个ID转换咯
回复

使用道具 举报

1

主题

12

帖子

128

积分

注册会员

Rank: 2

积分
128
发表于 2017-3-28 11:34:00 | 显示全部楼层
wly0410 发表于 2017-3-28 08:57
methylKit可以读入bismark的输出SAM文件,但是需要一定处理就是排序,这里借助SAMtools就可以了。samtools  ...

嗯嗯,我使用了这个了,已经把sam文件直接读进来了,我想知道能不能直接读extract的结果,来分析?
回复 支持 1 反对 0

使用道具 举报

1

主题

10

帖子

257

积分

中级会员

Rank: 3Rank: 3

积分
257
 楼主| 发表于 2017-2-24 00:40:25 | 显示全部楼层
3.2 bismark --bowtie2 -p 30 --score_min L,0,-0.6 --phred33-quals -N 1 /home/wenhao/BS-seq-temperature-muscle -1 D2-TE22_ATCACG_L008_R1.clean.trim.fastq.gz -2 D2-TE22_ATCACG_L008_R2.clean.trim.fastq.gz

应改成

bismark --bowtie2 -p 30 --score_min L,0,-0.6 --phred33-quals -N 1 /home/data/genome_file_dir -1 read1.fastq.gz -2 read2.fastq.gz

回复 支持 反对

使用道具 举报

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2017-2-24 09:08:53 | 显示全部楼层
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

10

主题

52

帖子

559

积分

版主

Rank: 7Rank: 7Rank: 7

积分
559
QQ
发表于 2017-2-28 23:03:44 | 显示全部楼层
很好的学习帖,可以参考坛主的链接使用插入代码,使代码更清晰
回复 支持 反对

使用道具 举报

1

主题

12

帖子

128

积分

注册会员

Rank: 2

积分
128
发表于 2017-3-23 11:51:09 | 显示全部楼层
楼主分析到那一步了?我现在已经完成map了,Bismark_methylation_extact也已经做完了,在研究怎么用methylKit分析后面的数据。一起交流一下,我现在遇到一个问题,怎么用methylKit读入,Bismark_extact出来的文件。你知道怎么操作嘛?
回复 支持 反对

使用道具 举报

0

主题

19

帖子

1166

积分

金牌会员

Rank: 6Rank: 6

积分
1166
发表于 2017-3-28 08:53:25 | 显示全部楼层
mehthykit网上有详细的教程说明
回复 支持 反对

使用道具 举报

0

主题

19

帖子

1166

积分

金牌会员

Rank: 6Rank: 6

积分
1166
发表于 2017-3-28 08:57:14 | 显示全部楼层
methylKit可以读入bismark的输出SAM文件,但是需要一定处理就是排序,这里借助SAMtools就可以了。samtools view -bS PE_1.fq_bismark_bt2_pe.sam >PE_1.bam
samtools sort PE_1.bam PE_1.sorted
samtools view -h PE_1.sorted.bam >PE_1.sorted.sam
回复 支持 反对

使用道具 举报

0

主题

19

帖子

1166

积分

金牌会员

Rank: 6Rank: 6

积分
1166
发表于 2017-3-28 18:42:14 | 显示全部楼层
这个需要格式转换才行
回复 支持 反对

使用道具 举报

0

主题

19

帖子

1166

积分

金牌会员

Rank: 6Rank: 6

积分
1166
发表于 2017-3-28 18:43:14 | 显示全部楼层
methylkit有自己一套分析流程,问什么要用extract数据?
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-10-20 07:25 , Processed in 0.036760 second(s), 31 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.