前言

最近一直在看小麥表觀遺傳相關的文獻，然后下載了文章中的數(shù)據(jù)進行分析，分析的流程大致跟網(wǎng)上教程相似，但是由于小麥基因組比較大（16G），研究不像人類，小鼠及一些模式植物廣泛，完善，所以中間有些分析內(nèi)容還是有些不同，于是把前期分析的代碼放在這里，以供交流學習
數(shù)據(jù)出處：https://doi.org/10.1186/s13059-020-01998-1

質(zhì)控比對

按照文章中的標準進行過濾和比對，因為樣品比較少，所以沒有批處理，也更直觀一些，中間可能挑了單個樣品跑了測試，所以有的樣品沒有出現(xiàn)在主腳本中

文章處理流程

#QC trim
java -jar /usrdata/users/hwwang/software/Trimmomatic-0.36/trimmomatic-0.36.jar PE -phred33 -threads 16 H3K27me3_1.fastq.gz H3K27me3_2.fastq.gz H3K27me3_1_trimmed.fastq.gz H3K27me3_1_unpaired.fastq.gz H3K27me3_2_trimmed.fastq.gz H3K27me3_2_unpaired.fastq.gz -validatePairs ILLUMINACLIP:/usrdata/users/hwwang/software/Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10 LEADING:5 TRAILING:5 MINLEN:20 2>>trim.log
java -jar /usrdata/users/hwwang/software/Trimmomatic-0.36/trimmomatic-0.36.jar PE -phred33 -threads 16 H3K36me3_1.fastq.gz H3K36me3_2.fastq.gz H3K36me3_1_trimmed.fastq.gz H3K36me3_1_unpaired.fastq.gz H3K36me3_2_trimmed.fastq.gz H3K36me3_2_unpaired.fastq.gz -validatePairs ILLUMINACLIP:/usrdata/users/hwwang/software/Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10 LEADING:5 TRAILING:5 MINLEN:20 2>>trim.log
java -jar /usrdata/users/hwwang/software/Trimmomatic-0.36/trimmomatic-0.36.jar PE -phred33 -threads 16 H3K4me3_1.fastq.gz H3K4me3_2.fastq.gz H3K4me3_1_trimmed.fastq.gz H3K4me3_1_unpaired.fastq.gz H3K4me3_2_trimmed.fastq.gz H3K4me3_2_unpaired.fastq.gz -validatePairs ILLUMINACLIP:/usrdata/users/hwwang/software/Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10 LEADING:5 TRAILING:5 MINLEN:20 2>>trim.log
java -jar /usrdata/users/hwwang/software/Trimmomatic-0.36/trimmomatic-0.36.jar PE -phred33 -threads 16 H3K9ac_1.fastq.gz H3K9ac_2.fastq.gz H3K9ac_1_trimmed.fastq.gz H3K9ac_1_unpaired.fastq.gz H3K9ac_2_trimmed.fastq.gz H3K9ac_2_unpaired.fastq.gz -validatePairs ILLUMINACLIP:/usrdata/users/hwwang/software/Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10 LEADING:5 TRAILING:5 MINLEN:20 2>>trim.log
java -jar /usrdata/users/hwwang/software/Trimmomatic-0.36/trimmomatic-0.36.jar PE -phred33 -threads 16 RNAPII_1.fastq.gz RNAPII_2.fastq.gz RNAPII_1_trimmed.fastq.gz RNAPII_1_unpaired.fastq.gz RNAPII_2_trimmed.fastq.gz RNAPII_2_unpaired.fastq.gz -validatePairs ILLUMINACLIP:/usrdata/users/hwwang/software/Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10 LEADING:5 TRAILING:5 MINLEN:20 2>>trim.log
#align
ref=/usrdata/users/hwwang/hychao/genome/iwgsc_refseqv1.0_all_chromosomes/bowtie2_index/wheat_bowtie_index
bowtie2 -p 12 --very-sensitive -x $ref -1 H3K27me3_1_trimmed.fastq.gz -2 H3K27me3_2_trimmed.fastq.gz > H3K27me3.sam
bowtie2 -p 12 --very-sensitive -x $ref -1 H3K36me3_1_trimmed.fastq.gz -2 H3K36me3_2_trimmed.fastq.gz > H3K36me3.sam
bowtie2 -p 12 --very-sensitive -x $ref -1 H3K4me3_1_trimmed.fastq.gz -2 H3K4me3_2_trimmed.fastq.gz > H3K4me3.sam
bowtie2 -p 12 --very-sensitive -x $ref -1 H3K9ac_1_trimmed.fastq.gz -2 H3K9ac_2_trimmed.fastq.gz > H3K9ac.sam
bowtie2 -p 12 --very-sensitive -x $ref -1 Input_1_trimmed.fastq.gz -2 Input_2_trimmed.fastq.gz > Input.sam
bowtie2 -p 12 --very-sensitive -x $ref -1 RNAPII_1_trimmed.fastq.gz -2 RNAPII_2_trimmed.fastq.gz > RNAPII.sam
#sort
samtools sort -@ 12 -o H3K27me3.bam H3K27me3.sam
samtools sort -@ 12 -o H3K36me3.bam H3K36me3.sam
samtools sort -@ 12 -o H3K4me3.bam H3K4me3.sam
samtools sort -@ 12 -o H3K9ac.bam H3K9ac.sam
samtools sort -@ 12 -o Input.bam Input.sam
samtools sort -@ 12 -o RNAPII.bam RNAPII.sam
#filter
sambamba view -h -t 8 -f bam -F 'not unmapped and not duplicate and mapping_quality >= 10' H3K27me3.bam > H3K27me3_aln.bam
sambamba view -h -t 8 -f bam -F 'not unmapped and not duplicate and mapping_quality >= 10' H3K36me3.bam > H3K36me3_aln.bam
sambamba view -h -t 8 -f bam -F 'not unmapped and not duplicate and mapping_quality >= 10' H3K4me3.bam > H3K4me3_aln.bam
sambamba view -h -t 8 -f bam -F 'not unmapped and not duplicate and mapping_quality >= 10' H3K9ac.bam > H3K9ac_aln.bam
sambamba view -h -t 8 -f bam -F 'not unmapped and not duplicate and mapping_quality >= 10' Input.bam > Input_aln.bam
sambamba view -h -t 8 -f bam -F 'not unmapped and not duplicate and mapping_quality >= 10' RNAPII.bam > RNAPII_aln.bam

生成比對結果

下面需要注意的一點就是samtools建立index生成 .bai 文件對染色體的長度是有限制的，因為小麥染色體長度在7，800M左右，運行的時候會報錯，所以在后面加一個-c參數(shù)，生成.csi文件

#qc_bam
ls  *aln.bam  |xargs -i samtools index -c {} 
ls  *aln.bam  | while read id ;do (nohup samtools flagstat $id > $(basename $id "aln.bam").stat & );done
grep N/A *.stat|grep %

IGV可視化

下一步就是生成IGV可視化的文件，這一步有一點小坑，就是bamCoverage軟件在進行歸一化的時候有一點改動，下面是之前的代碼

#igv calculate reads density
bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeTo1x 14547261565 -b H3K27me3_aln.bam -o H3K27me3.bigwig
bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeTo1x 14547261565 -b H3K36me3_aln.bam -o H3K36me3.bigwig
bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeTo1x 14547261565 -b H3K4me3_aln.bam -o H3K4me3.bigwig
bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeTo1x 14547261565 -b H3K9ac_aln.bam -o H3K9ac.bigwig
bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeTo1x 14547261565 -b Input_aln.bam -o Input.bigwig
bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeTo1x 14547261565 -b RNAPII_aln.bam -o RNAPII.bigwig

這里是我修改的

bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeUsing RPGC --effectiveGenomeSize 14547261565 -b H3K36me3_aln.bam -o H3K36me3.bigwig 2>>run.log
bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeUsing RPGC --effectiveGenomeSize 14547261565 -b H3K4me3_aln.bam -o H3K4me3.bigwig 2>>run.log
bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeUsing RPGC --effectiveGenomeSize 14547261565 -b H3K9ac_aln.bam -o H3K9ac.bigwig 2>>run.log
bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeUsing RPGC --effectiveGenomeSize 14547261565 -b Input_aln.bam -o Input.bigwig 2>>run.log
bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeUsing RPGC --effectiveGenomeSize 14547261565 -b RNAPII_aln.bam -o RNAPII.bigwig 2>>run.log

總的來說，跟文章中的處理過程是一致的

文章處理流程

最后就是call peak了！

#call peaks(histone)
macs2 callpeak -f BAM --nomodel --to-large -p 0.01 -broad -g 17e9 --bw 300 -c Input_aln.bam -t -n H3K27me3_result --outdir  ./
macs2 callpeak -f BAM --nomodel --to-large -p 0.01 -broad -g 17e9 --bw 300 -c Input_aln.bam -t -n H3K36me3_result --outdir  ./
macs2 callpeak -f BAM --nomodel --to-large -p 0.01 -broad -g 17e9 --bw 300 -c Input_aln.bam -t -n H3K4me3_result --outdir  ./
macs2 callpeak -f BAM --nomodel --to-large -p 0.01 -broad -g 17e9 --bw 300 -c Input_aln.bam -t -n H3K9ac_result --outdir  ./
#call peaks RNAPII
macs2 callpeak -f BAMPE --nomodel –q 0.001 --broad-cutoff 0.01 -g 17e9 --bw 300 -c Input_aln.bam -t -n RNAPIIRNAPII_result --outdir ./

上面的call peak代碼根據(jù)文章所述

文章分析流程

后面的分析內(nèi)容先簡單放點吧，因為還沒有全部完成，挑了一個數(shù)據(jù)，做了peak注釋

setwd("E:\\R_file\\test_chip\\")
library("ChIPseeker")
library("GenomicFeatures")
wheat_txdb <- loadDb("E:\\R_file\\wheat.sqlite")
RNAPII <- readPeakFile("H3K27me3_result_peaks.broadPeak.fold3")
peakAnno <- annotatePeak(RNAPII,tssRegion=c(-3000, 3000),TxDb=wheat_txdb)
plotAnnoPie(peakAnno)
plotAnnoBar(peakAnno)
vennpie(peakAnno)
upsetplot(peakAnno)

隨便幾個圖

更新一下

輸出注釋的peaks，一些圖可以嘗試畫一下

png("covplot_H3K27me3.png")
covplot(H3K27me3,weightCol=5)
dev.off()
##get promoter region +-3Kb and write out
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
write.table(promoter,"wheat_promoter.txt",row.names=FALSE,sep="\t")
#tagMatrix <- getTagMatrix(H3K27me3, windows=promoter)
#png("tagHeatmap.png")
#tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")
#dev.off()

##peak anno
peakAnno <- annotatePeak(H3K27me3,tssRegion=c(-3000, 3000),TxDb=txdb)
anno_out=as.data.frame(peakAnno)
write.table(anno_out,"H3K27me3_peakAnno.txt",row.names=FALSE,sep="\t")

標記數(shù)目比較多的話，可以使用 chromHMM 軟件分析不同標記組合的染色質(zhì)狀態(tài)

10種states

distribution

寫在最后

總的來說，前面的分析并不麻煩，就是有一些地方容易踩坑，寫下這些，主要是為了方便回看學習，也為了后面有做這方面的同學不要再犯跟我一樣的錯誤。關于小麥的 txDb 文件，怎樣制作網(wǎng)上都有教程，搜一下就可以了