一个MeDIP-seq实战（优秀学徒成果展）

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>

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

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

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

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 &