参考文献：The landscape of accessible chromatin in mammalian preimplantation embryos. Nature 2016 Jun 30;534(7609):652-7. PMID: 27309802
参考教程：https://www.jianshu.com/p/5bce43a537fd
https://www.jianshu.com/p/09e05bcd6981
miniconda
# https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/
wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh
source ~/.bashrc
#设置镜像。
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes
#other commands:https://www.jianshu.com/p/a20bb8c0bc5c
conda create -n atac -y python=2
conda info --envs
source activate atac
# 保证所有的软件都是安装在 atac这个环境下面
conda install -y sra-tools
conda install -y trim-galore samtools bedtools
conda install -y deeptools homer meme
conda install -y macs2 bowtie bowtie2
conda install -y multiqc
conda install -y sambamba
download sra
#config.sra
2-cell-1 SRR2927015
2-cell-2 SRR2927016
2-cell-4 SRR2927018
#prefetch
mkdir {sra,raw,clean,align,peaks,motif,qc}
cd ../sra
source activate atac
cat config.sra |while read id;do
arr=($id)
srr=${arr[1]}
nohup prefetch $srr &
done
## 默认下载目录：~/ncbi/public/sra/
fastq-dump
cat config.sra |while read id;
do arr=($id); srr=${arr[1]}; sample=${arr[0]}
nohup fastq-dump -A $sample -O ../raw --gzip --split-e ../$srr.sra &
done
QC
#config.raw
ls ../*_1.fastq.gz > 1 #complete path
ls ../*_2.fastq.gz > 2
paste 1 2 > config.raw
mkdir -p clean
cat config.raw |while read id;
do echo $id
arr=($id)
fq2=${arr[2]}
fq1=${arr[1]}
nohup trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paired -o clean $fq1 $fq2 &
done
cd ../qc
mkdir -p clean
fastqc -t 5 ../clean/2-cell-*gz -o clean
mkdir -p raw
fastqc -t 5 ../raw/2-cell-*gz -o raw
mapping
时间有限，只做了前两个数据
#config.clean
ls ../_1*gz >1
ls ../_2*gz >2
ls ../_1*gz | cut -d"/" -f 6 | cut -d"_" -f 1 >0
#/home/zhangjianing/nature_2016/clean/2-cell-1_1_val_1.fq.gz
paste 0 1 2 >config.clean
cat config.clean |while read id;
do echo $id; arr=($id); fq2=${arr[2]}; fq1=${arr[1]}; sample=${arr[0]}
bowtie2 -p 10 --very-sensitive -X 2000 -x ~/nature_2016/ref/mm10 -1 $fq1 -2 $fq2 |samtools sort -O bam -@ 5 -o - > ${sample}.raw.bam
samtools index ${sample}.raw.bam
bedtools bamtobed -i ${sample}.raw.bam > ${sample}.raw.bed
samtools flagstat ${sample}.raw.bam > ${sample}.raw.stat
sambamba markdup --overflow-list-size 600000 --tmpdir='./' -r ${sample}.raw.bam ${sample}.rmdup.bam
samtools index ${sample}.rmdup.bam
samtools flagstat ${sample}.rmdup.bam > ${sample}.rmdup.stat
##Calculate %mtDNA:( ref:https://www.biostars.org/p/170294/ )
mtReads=$(samtools idxstats ${sample}.rmdup.bam | grep 'chrM' | cut -f 3);
totalReads=$(samtools idxstats ${sample}.rmdup.bam | awk '{SUM += $3} END {print SUM}')
echo '==>' ${sample} 'mtDNA Content:' $(bc <<< "scale=2;100*$mtReads/$totalReads")'%'
samtools view -h -f 2 -q 30 ${sample}.rmdup.bam |grep -v chrM |samtools sort -O bam -@ 5 -o - > ${sample}.last.bam
samtools index ${sample}.last.bam
samtools flagstat ${sample}.la
st.bam > ${sample}.last.stat
bedtools bamtobed -i ${sample}.last.bam > ${sample}.bed
done
peak calling
ls *.bed | while read id ;
do macs2 callpeak -t $id -g mm --nomodel --shift -100 --extsize 200 -n ${id%%.*} --outdir ../peaks/ ;
# ${file%%.*}：拿掉第一个 “. ” 及其右边的字符串 (https://blog.csdn.net/chinalinuxzend/article/details/1826623)
done
计算插入片段长度
#config.fraction
last.bam path/2-cell-1 2-cell-1
last.bam path/2-cell-2 2-cell-2
#insert_fraction.R
cmd=commandArgs(trailingOnly=TRUE);
input=cmd[1]; output=cmd[2];
a=abs(as.numeric(read.table(input)[,1]));
png(file=output);
hist(a,
main="Insertion Size distribution",
ylab="Read Count",xlab="Insert Size",
xaxt="n",
breaks=seq(0,max(a),by=10)
);
axis(side=1,
at=seq(0,max(a),by=100),
labels=seq(0,max(a),by=100)
);
dev.off()
cat config.fraction | while read id;
do arr=($id); input=${arr[0]}; output=${arr[1]}
samtools view $input.last.bam | cut -f 9 > $output.txt #提取bam文件中的第9列（插入片段长度）
Rscript insert_fraction.R $output.txt $output.png;
done
2-cell-1
2-cell-2
FRIP value
bedtools intersect -a ../new/2-cell-1.bed -b 2-cell-1_peaks.narrowPeak |wc -l
148928
wc ../new/2-cell-1.bed
5105844
wc ../new/2-cell-1.raw.bed
5105844
ls *narrowPeak|while read id;
do
echo $id
bed=../new/$(basename $id "_peaks.narrowPeak").raw.bed
#ls -lh $bed
Reads=$(bedtools intersect -a $bed -b $id |wc -l|awk '{print $1}')
totalReads=$(wc -l $bed|awk '{print $1}')
echo $Reads $totalReads
echo '==> FRiP value:' $(bc <<< "scale=2;100*$Reads/$totalReads")'%'
done
overlap
#R studio serve安装失败，改用本地R studio
library(ChIPseeker)
library(ChIPpeakAnno)
setwd("/overlap")
list.files('./',"*.narrowPeak")
tmp = lapply(list.files('./',"*.narrowPeak"),function(x){
return(readPeakFile(file.path('./', x)))
})
`2-cell-` <- tmp #更改venn图中的名字，效果不是很好
ol <- findOverlapsOfPeaks(`2-cell-`[[1]],`2-cell-`[[2]])
png('overlapVenn.png')
makeVennDiagram(ol)
dev.off()
overlapVenn
6.deeptools可视化
cd /align
ls *last.bam |while read id;do
nohup bamCoverage -p 5 --normalizeUsing CPM -b $id -o ${id%%.*}.last.bw &
done
#TSS附近信号强度
cd /bw
ls *.bw | while read id;do
computeMatrix reference-point --referencePoint TSS -p 15 \
-b 10000 -a 10000 \
-R /ref/mm10.bed \
-S $id \
--skipZeros -o /TSS/${id%%.*}_TSS.gz \
--outFileSortedRegions ${id%%.*}_sortedRegions.bed;
#plot
plotHeatmap -m TSS/${id%%.*}_TSS.gz -out ${id%%.*}_TSS_Heatmap.png
plotProfile -m TSS/${id%%.*}_TSS.gz -out ${id%%.*}_TSS_Profile.png
(plotHeatmap -m ${id%%.*}_TSS.gz -out ${id%%.*}_TSS_Heatmap.pdf --plotFileFormat pdf --dpi 720
plotProfile -m ${id%%.*}_TSS.gz -out ${id%%.*}_TSS_Profile.pdf --plotFileFormat pdf --perGroup --dpi 720 )
#gene body 信号强度
ls *.bw | while read id;do
computeMatrix scale-regions -p 15 \
-R /ref/mm10.bed \
-S $id \
-b 10000 -a 10000 \
--skipZeros -o ${id%%.*}_body.gz
plotHeatmap -m ${id%%.*}_body.gz -out ${id%%.*}_body_Heatmap.png
plotProfile -m ${id%%.*}_body.gz -out ${id%%.*}_body_Profile.png
2-cell-2_TSS_Heatmap.png
2-cell-1_TSS_Heatmap.png
7.peaks注释
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("TxDb.Mmusculus.UCSC.mm10.knownGene")
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("org.Mm.eg.db")
library('CHIPpeakAnno')
library('TxDb.Mmusculus.UCSC.mm10.knownGene')
library('clusterProfiler')
setwd("/")
peakfiles = lapply(list.files('./',"*.bed"),function(x){return(readPeakFile(file.path('./', x)))})
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
peakAnno <- annotatePeak(peakfiles[[1]],tssRegion=c(-3000, 3000),TxDb=txdb, annoDb="org.Mm.eg.db")
peakAnno_df <- as.data.frame(peakAnno)
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
tagMatrix <- getTagMatrix(peakfiles[[1]], windows=promoter)
# 查看这些peaks在所有基因的启动子附近的分布情况，热图模式
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")
# 查看这些peaks在所有基因的启动子附近的分布情况，信号强度曲线图
plotAvgProf(tagMatrix, xlim=c(-3000, 3000), xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency")
2-cell-1
#others
plotAnnoPie(peakAnno)
plotAnnoBar(peakAnno)
#peak离最近基因的距离分布
plotDistToTSS(peakAnno, title="Distribution of transcription factor-binding loci\nrelative to TSS")
图片.png
图片.png
图片.png
#多图展示
peakAnnoList <- lapply(peakfiles, annotatePeak,TxDb=txdb,tssRegion=c(-3000, 3000))
plotAnnoBar(peakAnnoList)
plotDistToTSS(peakAnnoList)
图片.png
图片.png
homer peak annotation
# perl ~/miniconda3/envs/atac/share/homer-4.9.1-5/configureHomer.pl -install mm10
source activate atac
cd ../peaks
ls *.narrowPeak |while read id;
do
echo $id
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' $id > ${id%%.*}.homer_peaks.tmp
annotatePeaks.pl ${id%%.*}.homer_peaks.tmp mm10 \
1>${id%%.*}.peakAnn.xls \
2>${id%%.*}.annLog.txt
done
motifs
ls ../*.tmp |while read id;
do findMotifsGenome.pl $id mm10 ${id%%.*}_motifDir -len 8,10,12
done
图片.png
网友评论