ATACseq数据处理-补充/修正

crewkickse

已于 2023-07-05 15:37:21 修改

阅读量291

点赞数

文章标签： linux 学习

于 2023-07-03 13:19:51 首次发布

本文链接：https://blog.csdn.net/crewkickse/article/details/131405651

版权

fastqer-dump获得双端fastq文件后

过滤长度小于35的reads

nohup trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paired -o clean/ 1.fastq.gz 2.fastq.gz &

bowtie2比对参考基因组mm9：

nohup bowtie2 -p 5 --very-sensitive -X 2000 -x /share/home/wuqian/job/atac_test/cleandata/index/mm9 -1 SRR8528251_1.fastq -2 SRR8528251_2.fastq > SRR8528251.sam 2>&1 & 

##按照位置sort
samtools view -bS test.51.bam | samtools sort  -O bam  -@ 5 -o - > test.raw.sorted.bam

查看bam文件

samtools view SRR8528251_rmchrM_rmdup.bam|less -S

查看数据过滤前后reads数目的变化以及数据过滤结果

samtools flagstat SRR8528255_rmchrM.bam > SRR8528255_rmchrM.stat

.bai文件即为index文件

masc2 callpeak 命令改为：

nohup ls *.output.tagAlign.gz | while read id; do (macs2 callpeak -t $id -f BED -n "${id%%.*}" -g mm --shift -100 --extsize 200 --nomodel) ;done &

插入片段统计及FRiP值计算

FRiP：fraction of reads in called peak regions ;计算需要使用去除PCR重复之前的bed文件

peaks内的reads才是有效reads，落在peaks外的reads是无效reads，因此需要计算落在peaks内的reads占所有sam文件的reads的比例们就是所谓的FRiP值。

##bam转bed文件
nohup less /share/home/wuqian/job/shared/SRR_Acc_List.txt |while read id; do bedtools bamtobed -i ${id}.bam >${id}.raw.bed done 2>&1 &

##插入片段统计
samtools view SRR8528251_rmchrM_rmdup.bam |cut -f 9 >51.txt
##回到R语言

使用IDR(Irreproducible Discovery Rate)对不同重复生物学样本之间的差异peaks 进行提取
参考：CHIP-seq流程学习笔记（9）-使用IDR 软件对生物学重复样本间的差异peak进行提取_idr peak_垚垚爸爱学习的博客-CSDN博客

idr -s SRR8528251_peaks.narrowPeak SRR8528255_peaks.narrowPeak
wc idrValues.txt ##查看生成结果
##  67406 1348120 9216033 idrValues.txt 第一个数据就是重复reads数
## Number of reported peaks - 67406/67406 (100.0%)
## Number of peaks passing IDR cutoff of 0.05 - 32355/67406 (48.0%)

FRiP:

#生成bed文件
nohup ls *.nodup.bam|while read id;do (bedtools bamtobed -i $id >${id%%.*}.nodup.bed) ;done &
#批量计算FRiP
nohup less /share/home/wuqian/job/atac_test/atac_test.txt | while read id; do
echo $id; bed=${id}.last.bed; 
Reads=$(bedtools intersect -a $bed -b /share/home/wuqian/job/atac_test/mapping_genome/callpeak/${id}_summits.bed |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 2>&1

1.如何下载参考序列：

得到一个物种所有基因的TSS(转录起始位点)区域的bed文件 - 简书

2.bam转bw（使用deeptools bamCoverage）：

deeptools应用：- BAM和bigWig文件处理 - 质量控制 - 热图和其他描述性作图 - 其他

3.该步骤不需要经过Tn5偏移之后的文件，因为只有macs2 callpeak才需要，所以可以使用last.bam文件转换为bw文件之后进行可视化处理。

4.使用R语言提取无核小体、单核小体、双核小体和三核小体区域，并使用deeptools绘制热图和折线图(好慢)。

for (i in 1:length(read_51)) {
  read <- read_51[i]
  fragment_size <- length(read)
  if (fragment_size < fragment_sizes[1]) {
    nucleosome_free51 <- c(nucleosome_free51, read)
  } else if (fragment_sizes[2] <= fragment_size && fragment_size <= fragment_sizes[3]) {
    mononucleosome51 <- c(mononucleosome51, read)
  } else if (fragment_sizes[4] < fragment_size && fragment_size <= fragment_sizes[5]) {
    dinucleosome51 <- c(dinucleosome51, read)
  } else if (fragment_sizes[6] < fragment_size && fragment_size <= fragment_sizes[7]) {
    trinucleosome51 <- c(trinucleosome51, read)
  }

5.narrowPeak文件就是bedPeakFile

CHIPseeker用于peaks基因注释：

##从gff文件构建基因组注释db文件
library(GenomicFeatures)
metadata <- data.frame(name="Resource URL", value="https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/635/GCF_000001635.27_GRCm39/GCF_000001635.27_GRCm39_genomic.gff.gz")
txdb<- makeTxDbFromGFF("/share/home/wuqian/job/atac_test/genome/GCF_000001635.27_GRCm39_genomic.gff",
                       format="gff", organism="Mus musculus", taxonomyId=9940, dataSource="ncbigenomes", metadata=metadata)
saveDb(txdb,file = "TxDb.Mus.ncbi.GRCm39_genomic.sqlite")
txdb<-loadDb("TxDb.Mus.ncbi.GRCm39_genomic.sqlite")
makeTxDbPackage(txdb, version="1.0", maintainer="qian wu <wq_00211@163.com>", author="qian wu", destDir=".",license="Artistic-2.0",
                pkgname="TxDb.Mus.ncbi.GRCm39.genomicv1")
##Linux 下 TxDb.Oaries.Ensembl.Rambouilletv1　同级目录运行命令：
##R CMD build ./TxDb.Oaries.Ensembl.Rambouilletv1 

##peaks的基因注释
library(ChIPseeker)
library(rtracklayer)
library(org.Mm.eg.db)
library(TxDb.Mus.ncbi.GRCm39.genomicv1)
txdb<-TxDb.Mus.ncbi.GRCm39.genomicv1
peak51<-readPeakFile("SRR8528251_peaks.narrowPeak")
peak55<-readPeakFile("SRR8528255_peaks.narrowPeak")
##提取我们想要的染色体,小鼠19+XY
keepChr<-grepl("NC_",seqlevels(peak51))
seqlevels(peak51,pruning.mode="coarse")<-seqlevels(peak51)[keepChr]

seqlevels(peak55,pruning.mode="coarse")<-seqlevels(peak55)[keepChr]
seqlevels(mm39.refseq.db)
mm39.refseq.db
annotate51<-annotatePeak(peak51,tssRegion = c(-1000,1000),TxDb = txdb,annoDb = "org.Mm.eg.db")
annotate55<-annotatePeak(peak55,tssRegion = c(-1000,1000),TxDb = txdb,annoDb = "org.Mm.eg.db")
##查看一下注释结果
library(dplyr)
as.GRanges(annotate51) %>% head(3)
annotatepeaks_51<-as.data.frame(annotate51)
annotatepeaks_55<-as.data.frame(annotate55)
plotAnnoPie(annotate51)##饼图
plotAnnoPie(annotate55)
##可视化多个样本peak在1kb范围内的分布热图
files <- list.files(pattern = "*.narrowPeak")
files_list <- as.list(files)
names(files_list) <- c("SRR8528251","SRR8528255")
peakHeatmap(files_list, weightCol="V5", TxDb=txdb, 
            upstream=1000, downstream=1000, 
            color=rainbow(length(files)))

##TSS附近信号强度曲线图
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
tagMatrixList <- lapply(files_list, getTagMatrix, 
                        windows=promoter)
plotAvgProf(tagMatrixList,xlim=c(-3000,3000),
             xlab="Genomic Region (5'->3')", 
             ylab = "Read Count Frequency")

motif富集分析：

##UCSC官网下载对应的基因组注释文件，一定要是从一开始处理fastq文件开始就对应的注释文件，这里错误地使用mm39代替GRCm39，导致浪费了很多时间
 
参数：findMotifsGenome.pl <peak/bed file> <genome file> <output directory>
 
nohup findMotifsGenome.pl /share/home/wuqian/job/atac_test/callpeak/SRR8528255_peaks.narrowPeak /share/home/wuqian/job/atac_test/genome/GCF_000001635.27_GRCm39_genomic.fna motif_enrichment -preparse > motif.log 2>&1 &

TF footprinting(HINT-ATAC):

更新一下，目前这个软件所包含的注释基因组已经包括了hg19、hg38、mm10、mm39、mm9、zv10、zv9，大部分教程文章里没有提及mm39，以至于我学习的时候吓死了，差点以为又要自己手动构建参考基因组文件（为什么是又。。。辛酸泪，以后比对之前一定要选好参考基因组文件，推荐Ensmembl）

##安装软件
pip install --user cython numpy scipy
pip install --user RGT --no-binary RGT

##进入home目录下的rgtdata文件夹
cd ~/rgtdata
python setupGenomicData.py --mm39

##TF footprinting
##输入文件需要包含一个bam文件(indexed and sorted)和callpeak之后生成的narrowpeak文件
nohup rgt-hint footprinting --atac-seq --paired-end --organism=mm39 --output-location=./ --output-prefix=55 ../mapping_genome/rmdup/SRR8528255_rmchrM_rmdup.bam ../callpeak/SRR8528255_peaks.narrowPeak 2>&1 &

这个软件只考虑peak区域进行footprinting，每个命令结束后将在当前文件下生成一个包含footprinting基因组位置的.bed输出文件，以及一个以.info结尾的文件，包含来自库的统计信息,即总提取次数等。

参考HINT-ATAC官网教程：Introduction — RGT documentation

使用HINT-ATAC进行ATAC-Seq的footprinting分析 - 简书