ATAC-Seq基础分析+高级分析+多组学分析

# 提取100w条reads代码
ls *gz|while read id;do zcat $id|head -4000000|gzip > ./100/${id};done

wget https://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_chromosome_files/TAIR10_chr_all.fas
wget https://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_gff3/TAIR10_GFF3_genes_transposons.gff
# 改名
mv TAIR10_chr_all.fas TAIR10_genome.fa
mv TAIR10_GFF3_genes_transposons.gff TAIR10.gff3

conda install -c bioconda -y fastqc
conda install -c bioconda -y trimmomatic
conda install -c bioconda -y bwa
conda install -c bioconda -y samtools
conda install -c bioconda -y picard
conda install -c bioconda -y macs2
conda install -c bioconda -y bedtools
conda install -c bioconda -y deeptools
conda install -c bioconda -y homer

fastqc -t 8 -o ./ *.fastq.gz

# 定义一些变量，后面修改起来方便
filename=SRR7512044
 
trimmomatic PE -threads 10 
./${filename}_1.fastq.gz ./${filename}_2.fastq.gz 
./${filename}_1.clean.fq.gz ./${filename}_1.drop.fq.gz 
./${filename}_2.clean.fq.gz ./${filename}_2.drop.fq.gz 
HEADCROP:15 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:24

bwa index -p /路径/bwaIndexTair /路径/TAIR10_genome.fa

# 定义一些变量，后面用
filename=SRR7512044
index=/路径/bwaIndexTair
 
bwa mem -M -t 8 
-R "@RGtID:${filename}tSM:${filename}tLB:WXStPL:Illumina" 
$index 
${filename}_1.clean.fq.gz 
${filename}_2.clean.fq.gz 
| samtools sort -O bam -@ 10 -o ./${filename}.raw.bam
#无论干啥，生成bam文件后就要建立bam索引，并检查比对情况
samtools index -@ 10 -b ./${filename}.raw.bam
samtools flagstat ./${filename}.raw.bam > ./${filename}.raw.stat
cat ${filename}.raw.stat

14608455   0 in total (QC-passed reads   QC-failed reads) reads总数
37967   0 secondary                                       出现比对到参考基因组多个位置的reads数
0   0 supplementary                                       可能存在嵌合的reads数
0   0 duplicates                                          重复的reads数
14590894   0 mapped (99.88% : N/A)                        比对到参考基因组上的reads数
14570488   0 paired in sequencing                         属于PE read的reads总数。
7285244   0 read1                                         PE read中Read 1 的reads 总数。
7285244   0 read2                                         PE read中Read 2 的reads 总数。
14507068   0 properly paired (99.56% : N/A)               完美比对的reads总数。PE两端reads比对到同一条序列，且根据比对结果推断的插入片段大小符合设置的阈值。
14551500   0 with itself and mate mapped                  PE两端reads都比对上参考序列的reads总数。
1427   0 singletons (0.01% : N/A)                         PE两端reads，其中一端比上，另一端没比上的reads总数。
26260   0 with mate mapped to a different chr             PE read中，两端分别比对到两条不同的序列的reads总数。
17346   0 with mate mapped to a different chr (mapQ>=5)   PE read中，两端分别比对到两条不同

picard MarkDuplicates 
REMOVE_DUPLICATES=true 
I=${filename}.raw.bam 
O=${filename}.rmdup.bam 
M=${filename}.rmdup.log
# 实时监测我们的数据发生了什么变化
samtools index ${filename}.rmdup.bam
samtools flagstat ${filename}.rmdup.bam > ./${filename}.rmdup.stat

samtools view -h -f 2 -q 30 ${filename}.rmdup.bam 
| grep -v -e "mitochondria" -e "*" -e "chloroplast" 
| samtools sort -O bam -@ 10 -o - > ${filename}.last.bam
# 实时监测我们的数据发生了什么变化
samtools index ${filename}.last.bam
samtools flagstat ${filename}.last.bam > ./${filename}.last.stat

picard CollectInsertSizeMetrics 
I=./${filename}.last.bam 
O=./${filename}.last.insertsize 
H=./${filename}.last.hist.pdf

effective_genome_size=119481543
 
bedtools bamtobed -i ${filename}.last.bam > ${filename}.last.bed
macs2 callpeak 
-t ${filename}.last.bed 
-g ${effective_genome_size} 
--nomodel --shift -100 --extsize 200 
-n ${filename}.last 
-q 0.01 --outdir ./peaks

前几行是callpeak时的命令。后面是：
1.染色体号
2.peak起始位点
3.peak结束位点
4.peak区域长度
5.peak的峰值位点（summit position）
6.peak 峰值的高度（pileup height at peak summit, -log10(pvalue) for the peak summit）
7.peak的富集倍数（相对于random Poisson distribution with local lambda）

narrowPeak文件是BED6 4格式，可以上传到IGV查看
1.染色体号
2.peak起始位点
3.peak结束位点
4.peak name（样本名 _peak_1的后缀之类的）
5.-10*log10qvalue
6.正负链
7.fold change
8.-log10pvalue
9.-log10qvalue
10.相对于峰的峰顶位置（relative summit position to peak start）

BED格式的文件，包含peak的summits（峰顶）位置
1.染色体号
2.峰的峰顶起始位点（注意是峰顶，不是峰，他只是一个点）
3.峰的峰顶起始位点（因此二者就差一个碱基的位置）
4.峰顶（样本名 _peak_1的后缀之类的）
5.-log10pvalue

Reads=$(bedtools intersect -a ${filename}.last.bed -b ${filename}_peaks.narrowPeak |wc -l|awk '{print $1}')
totalReads=$(wc -l ${filename}.last.bed|awk '{print $1}') 
echo ${Reads} ${totalReads} ${filename} 'FRiP_score:' $(echo "scale=2;100*${Reads}/${totalReads}"|bc)'%'

bamCoverage --bam ${filename}.last.bam -o ${filename}.last.bw 
    --binSize 10 
    --normalizeUsing RPGC 
    --effectiveGenomeSize ${effective_genome_size} 
    --ignoreForNormalization chrX 
    --extendReads

computeMatrix reference-point --referencePoint TSS  -p 15 
    -b 2000 -a 2000 
    -R ./refseq.bed 
    -S *.bw 
    -o TSS.gz 
    --skipZeros 
    --outFileSortedRegions Heatmap1sortedRegions.bed
plotHeatmap -m TSS.gz -out TSS.Heatmap.pdf --plotFileFormat pdf
plotProfile -m TSS.gz -out TSS.Profile.pdf --plotFileFormat pdf --perGroup

computeMatrix scale-regions  -p 15  
    -b 2000 -a 2000 
    -R ./refseq.bed 
    -S *.bw 
    --skipZeros -o body.gz
plotHeatmap -m body.gz -out body.Heatmap.pdf --plotFileFormat pdf
plotProfile -m body.gz -out body.Profile.pdf --plotFileFormat pdf --perGroup

multiBigwigSummary bins -b *.last.bw -o number.of.bins.npz

plotCorrelation -in number.of.bins.npz 
--corMethod pearson 
--skipZeros 
--plotTitle "Pearson Correlation of Average Scores Per Transcript" 
--plotFileFormat pdf 
--whatToPlot heatmap --colorMap RdYlBu --plotNumbers 
-o heatmap_PearsonCorr_bigwigScores.pdf 
--outFileCorMatrix PearsonCorr_bigwigScores.tab

plotCorrelation -in number.of.bins.npz 
--corMethod pearson 
--skipZeros 
--plotTitle "Pearson Correlation of Average Scores Per Transcript" 
--plotFileFormat pdf 
--whatToPlot scatterplot 
-o scatterplot_PearsonCorr_bigwigScores.pdf 
--outFileCorMatrix PearsonCorr_bigwigScores.tab

#if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install("clusterProfiler")
 
library(clusterProfiler)
library(ChIPseeker)
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("public/bioinfoYu/genome/tair10/TAIR10.gff",
                        format="gff")    #也可以使用gtf
keytypes(txdb)    #感兴趣的话，可以用以下方法探索txdb都包含了什么内容
keys(txdb)

#读入单个summits文件
peaks <- readPeakFile("SRR7512044_summits.bed")
#结构注释
peakAnno <- annotatePeak(peaks,
                         TxDb=txdb,
                         tssRegion=c(-1000, 1000))

#注释完，进行可视化，多种图可供选择
plotAnnoBar(peakAnno)
plotDistToTSS(peakAnno)
vennpie(peakAnno)
plotAnnoPie(peakAnno)
#install.packages("ggupset")
library(ggupset)
upsetplot(peakAnno)
#install.packages("ggimage")
library(ggimage)
upsetplot(peakAnno, vennpie=TRUE)
 
#最后将我们的注释结果转为数据框，便于查看
df <- as.data.frame(peakAnno)
#将注释到的基因提取出来（第14列），用于后续功能分析
gene <- df[,14]

#一次也可以读入多个summits文件，使用list存储，然后使用lapply注释
files = list(peak1 = "peaks.bed", peak2 = ("peak.bed"))
peakAnnoList <- lapply(files, 
                       annotatePeak,
                       TxDb=txdb,
                       tssRegion=c(-2000, 2000))
plotAnnoBar(peakAnnoList)
plotDistToTSS(peakAnnoList)

wget http://www.bioconductor.org/packages/release/data/annotation/src/contrib/org.At.tair.db_3.13.0.tar.gz
tar -zvxf org.At.tair.db_3.13.0.tar.gz

#从本地下载OrgDb库（如果你直接下载到了R的lib里面，就直接library也可）
install.packages("/public/bioinfoYu/genome/tair10/org.At.tair.db", repos=NULL)
#另外一种方法
#if (!requireNamespace("BiocManager", quietly = TRUE))install.packages("BiocManager")
#BiocManager::install("org.At.tair.db")
#载入库
library(org.At.tair.db)
#感兴趣的话，看看org.At.tair.db包含的都是些啥
#查看包含的列名
keytypes(org.At.tair.db)
#查看包含的基因ID
keys(org.At.tair.db)
#随便选择一个基因和一些列名，看一看到第是什么东东，不过因为导入了其他包，下面命令报错，我至今还没有解决
select(org.At.tair.db,keys="AT1G01010",columns = c("ENTREZID","GO","SYMBOL"))

#首先，去除版本号，如AT1G50040.1后面的.1，当然我们的测试集没有，因此不用做
#gene_rm <- gsub(".[0-9] $","",gene)
#下面进行ID转换
keytypes(org.At.tair.db)
gene_last <- bitr(geneID = gene, 
                    fromType = "TAIR",
                    toType = c("ENTREZID", "SYMBOL", "GENENAME"),
                    OrgDb = org.At.tair.db)
 
# geneID    需要转换的ID
# fromType  当前ID类型
# toType    转换成什么ID，使用keytypes()查看有哪些类型
# OrgDb     注释数据库

# 因为keyType不接受数据框类型的数据，因此我们转换为字符型
genelist <- gene.symbol[,2]  #注意我这里把ENTREZID单独拿出来了，所以下面自然是keyType = "ENTREZID",
go_result <- enrichGO(genelist,
                      keyType = "ENTREZID", 
                      OrgDb = org.At.tair.db, 
                      ont = "BP", #BP,MF,CC,ALL可选
                      pAdjustMethod = "BH", 
                      qvalueCutoff = 0.05, 
                      readable = FALSE)
#可视化，showCategory=20调整显示的条目多少
dotplot(go_result, showCategory=20)
#将结果转为数据框，可以选择导出
goresult <- as.data.frame(go_result)

#geneClusters必须是数据框类型，ID必须是ENTREZID类型
class(gene.symbol)
compKEGG <- compareCluster(geneClusters=gene.symbol, 
                           fun = "enrichKEGG",
                           organism = "ath",    #在 http://www.genome.jp/kegg/catalog/org_list.html 上查看物种简写
                           pvalueCutoff  = 0.05, 
                           pAdjustMethod = "BH")
dotplot(compKEGG, showCategory = 20, title = "KEGG Pathway Enrichment Analysis")

# 提取peaks的位置信息文件
awk '{print $4"t"$1"t"$2"t"$3"t "}' 
SRR7512044_peaks.narrowPeak > SRR7512044.homer_peaks.tmp
# 寻找motif
findMotifsGenome.pl SRR7512044.homer_peaks.tmp 
./TAIR10_genome.fa   #指定参考基因组
./SRR7512044.motif   #输出文件夹的名字
-size 200 -len 8,10,12

# 安装包
if (!requireNamespace("BiocManager", quietly = TRUE))install.packages("BiocManager")
BiocManager::install("DiffBind")
library(DiffBind)
#导入数据，dba函数会将会将bam文件一同导入，因此csv中的路径一定要准确
tamoxifen <- dba(sampleSheet="atDiff.csv")
#查看样本间的相关性
plot(tamoxifen)

#计算每个样本中的reads数
tamoxifen <- dba.count(tamoxifen)
#简单查看计算结果
info <- dba.show(tamoxifen)
libsizes <- cbind(LibReads=info$Reads, FRiP=info$FRiP, PeakReads=round(info$Reads * info$FRiP))
rownames(libsizes) <- info$ID
libsizes

#统计结果
       LibReads FRiP PeakReads
cmt2_1  6720462 0.29   1948934
cmt2_2  7747491 0.33   2556672
wild_1  7078309 0.36   2548191
wild_2  8030752 0.34   2730456
ddcc_1  6370232 0.30   1911069
ddcc_2  6956429 0.35   2434750

#标椎化结果
       FullLibSize  NormFacs NormLibSize
cmt2_1     6720462 0.9398443     7150612
cmt2_2     7747491 1.0834724     7150612
wild_1     7078309 0.9898885     7150612
wild_2     8030752 1.1230859     7150612
ddcc_1     6370232 0.8908652     7150612
ddcc_2     6956429 0.9728438     7150612

#分组，格式是表头在最前面，要分的组依次写在后面，只能两两比较，因此后面只能写两组，但可以多执行几次，都会追加到tamoxifen 中
#分组1，后面使用contrast=1单独查看
tamoxifen <- dba.contrast(tamoxifen, contrast=c("Factor","ddcc","wild"))
#分组2，后面使用contrast=2单独查看
tamoxifen <- dba.contrast(tamoxifen, contrast=c("Factor","cmt2","wild"))
#当然还可有分组3,4,5等，均可以使用contrast=number单独查看
tamoxifen
#按照分组分别进行差异分析，默认使用DESeq2进行计算，可以选择method = DBA_EDGER(edgR)，或者两个都要method = DBA_ALL_METHODS
tamoxifen <- dba.analyze(tamoxifen)
dba.show(tamoxifen, bContrasts=TRUE)
#查看差异分析的结果与导出为csv文件
#查看第1组差异分析的结果
tamoxifen.DB_1 <- dba.report(tamoxifen,contrast=1)
tamoxifen.DB_1
write.csv(tamoxifen.DB_1, file="diffBind_result.csv")
 
#查看第2组差异分析的结果
tamoxifen.DB_2 <- dba.report(tamoxifen,contrast=2) tamoxifen.DB_2 write.csv(tamoxifen.DB_2, file="diffBind_result.csv") ##传统的上下调在找差异peak中称为“Gain”或“Loss” #查看第1组的差异peak数量“Fold>0或<0”控制是“Gain”或“Loss” sum(tamoxifen.DB$Fold>0,contrast=1)
sum(tamoxifen.DB$Fold<0,contrast=1) #查看第2组的差异peak数量“Fold>0或<0”控制是“Gain”或“Loss” sum(tamoxifen.DB$Fold>0,contrast=2)
sum(tamoxifen.DB$Fold<0,contrast=2)

#韦恩图可视化
dba.plotVenn(tamoxifen, contrast=1, bDB=TRUE, bGain=TRUE, bLoss=TRUE, bAll=FALSE)
 
#PCA
dba.plotPCA(tamoxifen) #这里可以使用所有样本进行PCA
dba.plotPCA(tamoxifen, contrast=1, label=DBA_FACTOR)#单独对分组1进行PCA
 
#MA plots
dba.plotMA(tamoxifen, contrast=1)
 
#Volcano plots
dba.plotVolcano(tamoxifen, contrast=1)
 
#Box plots
dba.plotBox(tamoxifen, contrast=1)
 
#Heatmap
hmap <- colorRampPalette(c("red", "black", "green"))(n = 13)
#对所有样本的所有的差异peak画热图
dba.plotHeatmap(tamoxifen,correlations=FALSE,scale="row",colScheme = hmap)
#单独对分组1中所有的差异peak画热图
dba.plotHeatmap(tamoxifen, contrast=1, correlations=FALSE, scale="row", colScheme = hmap)
dba.plotHeatmap(tamoxifen,correlations=FALSE,scale="row",colScheme = hmap)

#Profiling and Profile Heatmaps
dba.plotProfile(tamoxifen)