其他
一个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文件
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
#注，这块命令目录和图片上的有出入，我后面移动了一下文件的位置
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 里面下载下面这个文件
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格式的图
导入电脑查看
若有多个样本，批处理脚本
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点评：
从生信工程师的角度来看，此实战经验分享可以说是满分了！
但是，作者毕竟是本科生，完全忽略了数据背后的生物学知识，所以并未重现文章的几个关键性图表及结论，这就是目前普遍的科研服务公司所以提供的标准流程服务与广大科研工作者需要的个性化分析之间的矛盾的最真实写照。
独家福利
如果需要组装自己的服务器；代办生物信息学服务器
如果需要帮忙下载海外数据(GEO/TCGA/GTEx等等)，点我？
如果需要线下辅导及培训，看招学徒
如果需要个人电脑：个人计算机推荐
如果需要置办生物信息学书籍，看：生信人必备书单
如果需要实习岗位：实习职位发布
如果需要售后：点我