一个MeDIP-seq实战(优秀学徒成果展)

前传: 一个简单进行SNP分析的实战例子(优秀学徒成果展)

一、安装软件

安装samtools(常用软件,可以安装在自己的主环境下)
cd ~/package
wget -c https://sourceforge.net/projects/samtools/files/samtools/1.8/samtools-1.8.tar.bz2
tar -xvf samtools-1.8.tar.bz2 -C ~/biosoft/
cd ~/biosoft/samtools-1.8
make
./samtools
ln -s ~/biosoft/samtools-1.8/samtools ~/biosoft/mybin/
下面新建一个conda环境来放处理MeDIP-seq所用到的软件
conda create -n MeDIP
source activate MeDIP
conda install macs2
conda install fastqc 
conda install bowtie2
conda install deeptools

二、下载数据

根据GSE号,查到原始数据仅有两个样本

方法一:可以选择直接单个下载
wget -c <ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP018/SRP018845/SRR764931/SRR764931.sra>

wget -c <ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP018/SRP018845/SRR764932/SRR764932.sra>
方法二:下载SraRunInfo.csv文件

img

cat SraRunInfo.csv | cut -f 10 -d ',' | grep SRR  > runinfo.csv
cat runinfo.csv | xargs -i echo wget -c {} \$ >>sra.sh
bash sra.sh 
数据转换
mkdir ~/MeDIP/fastq
ls *sra |while read id; do fastq-dump --split-3 $id -O ~/MeDIP-seq/fastq ;done
##注,这块命令目录和图片上的有出入,我后面移动了一下文件的位置

img

ls *fastq |xargs fastqc -t 8 -o ~/MeDIP-seq/qc

数据质量还是很好的

三、找peaks(四步曲)

3.1 比对

下载参考基因组的索引和注释文件

mkdir reference && cd reference
wget -4 -q ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip ##-4 使用IPv4协议
unzip mm10.zip
单个比对
bowtie2 -p 6 -x ~/MeDIP-seq/reference/mm10 -U SRR764931.fastq | samtools sort -O bam -o SRR764931.bam
多个样本比对
vim align.sh

bin_bowtie2=bowtie2
bowtie2_index=/home/yjzhang/MeDIP-seq/reference/mm10
ls ~/MeDIP-seq/fastq/*.fastq |while read id;
do 
file=$(basename $id )
sample=${file%%.*}
echo $file $sample
$bin_bowtie2  -p 8  -x  $bowtie2_index -U  $id | samtools sort  -O bam  -@ 8 -o - > ${sample}.bam 
done 

bash align.sh
比对所用脚本
统计比对率,结果挺高的
cd ~/MeDIP-seq/align
ls  *.bam  |xargs -i samtools index {} 
ls  *.bam  | while read id ;do (nohup samtools flagstat $id > ~/MeDIP-seq/qc/align-qc/$(basename $id ".bam").stat & );done

3.2 找peaks

对bam文件建索引
ls  *.bam  |xargs -i samtools index {} 
使用macs2进行找peaks
macs2 callpeak -t SRR764931.bam -m 10 30 -p 1e-5 -f BAM -g mm -n /home/yjzhang/MeDIP-seq/peak/SRR764931 2>SRR764931.masc2.log
报错,查看报错信息,缺少libgfortran.so.1 ,按理说用conda安装的macs2应该不会缺少配置的,查看安装macs2时的信息,只安装了libgfortran.so.3,问题可能出现在这里
尝试一下安装libgfortran.so.1
conda install libgfortran=1
再次执行macs命令,成功,输出4个结果文件,说明第一次执行不成功是libgfortran版本的问题
若有多个样本,可以用下面的脚本
cd ~/MeDIP-seq/align
vim macs.sh

source activate  MeDIP
ls  *.bam |cut -d"." -f 1 |while read id;
do 
    if [ ! -s /home/yjzhang/MeDIP-seq/peak/${id}_summits.bed ];
    then 
echo $id 
nohup macs2 callpeak -t $id.bam -m 10 30 -p 1e-5 -f BAM -g mm -n /home/yjzhang/MeDIP-seq/peak/$id 2> $id.log &  
    fi 
done  
source deactivate
运行bash macs.sh 结果如下

3.3 bam格式转换为bw格式

deeptools提供bamCoverage和bamCompare进行格式转换
bamCoverage -p 10 --normalizeUsing RPKM -b SRR764931.bam -o SRR764931.bw 
若有多个样本,批量转换
cd  ~/MeDIP-seq/align
vim bw.sh
source activate MeDIP
ls *.bam |while read id;do
nohup bamCoverage --normalizeUsing RPKM -p 10 -b $id -o ./${id%%.*}.bw & 
done 
source deactivate

3.4 查看TSS附件信号强度并作图

在执行这一步之前,要先下载小鼠的mm10.ucsc.refseq.bed文件
可以在UCSC的table browser 里面下载下面这个文件

img

mkdir TSS && cd TSS
##搜峰
computeMatrix reference-point --referencePoint TSS -b 10000 -a 10000 -R ~/MeDIP-seq/reference/mm10.ucsc.refseq.bed -S ~/MeDIP-seq/align/SRR764931.bw --skipZeros -o matrix1_SRR764931_TSS.gz
##画图
plotHeatmap -m matrix1_SRR764931_TSS.gz -out SRR764931.png
plotHeatmap -m matrix1_SRR764931_TSS.gz -out SRR764931.pdf --plotFileFormat pdf  --dpi 720
报错,numpy模块出错,查询了一下报错,貌似是numpy版本过低,升级一下版本
用pip install numpy –upgrade来更新,不行,因为还有用conda安装的别的软件也依赖numpy。
尝试一下用conda更新
conda install numpy (conda默认安装适合当前环境最新的版本)或者conda update numpy
版本从numpy-1.9.2改成numpy-1.14.5
接下来再次执行命令就能成功运行,而且我试了下也不影响其他软件运行
画png格式的图
画pdf格式的图
导入电脑查看

img

若有多个样本,批处理脚本
cd ~/MeDIP-seq/TSS/

vim tss.sh

bed=~/MeDIP-seq/reference/mm10.ucsc.refseq.bed
for id in  /home/yjzhang/MeDIP-seq/align/*.bw ;
do 
echo $id
file=$(basename $id )
sample=${file%%.*} 
echo $sample  

computeMatrix reference-point  --referencePoint TSS  -p 5  \
-b 10000 -a 10000    \
-R  $bed \
-S $id  \
--skipZeros  -o ./matrix1_${sample}_TSS.gz  \
## both plotHeatmap and plotProfile will use the output from     computeMatrix
plotHeatmap -m ./matrix1_${sample}_TSS.gz  -out ./${sample}_Heatmap.png
plotHeatmap -m ./matrix1_${sample}_TSS.gz  -out ./${sample}_Heatmap.pdf --plotFileFormat pdf  --dpi 720  
done 

## 使用命令提交脚本:
nohup bash tss.sh 1>10k.log &

参考: http://www.bio-info-trainee.com/2352.html http://www.bio-info-trainee.com/2494.html

jimmy点评:

从生信工程师的角度来看,此实战经验分享可以说是满分了!

但是,作者毕竟是本科生,完全忽略了数据背后的生物学知识,所以并未重现文章的几个关键性图表及结论,这就是目前普遍的科研服务公司所以提供的标准流程服务与广大科研工作者需要的个性化分析之间的矛盾的最真实写照。

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2018-10-07

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏小白课代表

Autodesk Revit 2014安装教程

Revit是Autodesk公司一套系列软件的名称。Revit系列软件是专为建筑信息模型(BIM)构建的,可帮助建筑设计师设计、建造和维护质量更好、能效更高的建...

50920
来自专栏IMWeb前端团队

[flash相关]crossBridge生成的库文件体积优化

本文作者:IMWeb 黄龙 原文出处:IMWeb社区 未经同意,禁止转载 不明白crossBridge是什么的可以看下这里 http://imweb....

23060
来自专栏安富莱嵌入式技术分享

【安富莱原创开源应用第3期】花式玩转网络摄像头之VNC远程桌面版本,稳定运行2年不死机

1、前段时间开源了一个网络摄像头的TCP版本 https://www.cnblogs.com/armfly/p/9173167.html,这次再来一个远程VNC...

14420
来自专栏用户2442861的专栏

使用事件驱动模型实现高效稳定的网络服务器程序

http://www.cnblogs.com/hnrainll/p/3625597.html

10410
来自专栏小白课代表

Autodesk Revit 2016安装教程

Revit是Autodesk公司一套系列软件的名称。Revit系列软件是专为建筑信息模型(BIM)构建的,可帮助建筑设计师设计、建造和维护质量更好、能效更高的建...

50440
来自专栏区块链源码分析

超级账本(Hyperledger Fabric)源码分析之一:总览

1)Go,注意设置好gopath(笔者安装的是go1.8.3,对应的源码是v1.0.0这个tag,版本不对可能会出现编译不过或者运行出现问题)

52650
来自专栏北京马哥教育

几种经典的网络服务器架构模型的分析与比较

前言 事件驱动为广大的程序员所熟悉,其最为人津津乐道的是在图形化界面编程中的应用;事实上,在网络编程中事件驱动也被广泛使用,并大规模部署在高连接数高吞吐量的服务...

29350
来自专栏用户2442861的专栏

Redis作者谈Redis应用场景

毫无疑问,Redis开创了一种新的数据存储思路,使用Redis,我们不用在面对功能单调的数据库时,把精力放在如何把大象放进冰箱这样的问题上,而是利用Redis...

30620
来自专栏ATYUN订阅号

腾讯开源围棋AI程序PhoenixGo,复现AlphaGo Zero

PhoenixGo是一个围棋AI程序,它执行AlphaGo Zero论文“掌握无人知识的Go游戏”。它也被称为FoxGo中的“BensonDarr”,CGOS中...

25420
来自专栏java一日一条

干货!如何正确使用Git Flow

我们已经从SVN 切换到Git很多年了,现在几乎所有的项目都在使用Github管理, 本篇文章讲一下为什么使用Git, 以及如何在团队中正确使用。

17240

扫码关注云+社区

领取腾讯云代金券