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 create -n MeDIP
source activate MeDIP
conda install macs2
conda install fastqc
conda install bowtie2
conda install deeptools
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>
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
数据质量还是很好的
下载参考基因组的索引和注释文件
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
ls *.bam |xargs -i samtools index {}
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
conda install libgfortran=1
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
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
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
conda install numpy (conda默认安装适合当前环境最新的版本)或者conda update numpy
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点评:
从生信工程师的角度来看,此实战经验分享可以说是满分了!
但是,作者毕竟是本科生,完全忽略了数据背后的生物学知识,所以并未重现文章的几个关键性图表及结论,这就是目前普遍的科研服务公司所以提供的标准流程服务与广大科研工作者需要的个性化分析之间的矛盾的最真实写照。