在R中也可以跑ATAC了，主要是上游分析也可以跑，惊讶把，来来，咱来跑一下，看一下咋样

1、安装需要的R包

这步是最麻烦的，大家耐心点，ChIPQC这个包好像有点问题，大家下的时候注意一下

## 使用西湖大学的 Bioconductor镜像
options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor")
options("repos"=c(CRAN="https://mirrors.westlake.edu.cn/CRAN/"))

# 本课程的包
install.packages('BiocManager')
BiocManager::install('RockefellerUniversity/RU_ATACseq',subdir='atacseq')

# CRAN and Bioconductor
install.packages('BiocManager')
BiocManager::install('methods')
BiocManager::install('ggplot2')
BiocManager::install('rmarkdown')
BiocManager::install('ShortRead')
BiocManager::install('ashr')
BiocManager::install('ChIPQC')
BiocManager::install('DiffBind')
BiocManager::install('BSgenome.Hsapiens.UCSC.hg19')
BiocManager::install('Rsubread')
BiocManager::install('Rbowtie2')
BiocManager::install('R.utils')
BiocManager::install('Rsamtools')
BiocManager::install('BSgenome.Hsapiens.UCSC.hg38')
BiocManager::install('rtracklayer')
BiocManager::install('ChIPseeker')
BiocManager::install('soGGi')
BiocManager::install('GenomicAlignments')
BiocManager::install('TxDb.Hsapiens.UCSC.hg19.knownGene')
BiocManager::install('DESeq2')
BiocManager::install('BSgenome.Mmusculus.UCSC.mm10')
BiocManager::install('TxDb.Hsapiens.UCSC.hg38.knownGene')
BiocManager::install('tracktables')
BiocManager::install('clusterProfiler')
BiocManager::install('TxDb.Mmusculus.UCSC.mm10.knownGene')
BiocManager::install('devtools')
BiocManager::install('tidyr')
BiocManager::install('DT')
BiocManager::install('dplyr')
BiocManager::install('rGREAT')
BiocManager::install('MotifDb')
BiocManager::install('Biostrings')
BiocManager::install('GenomicRanges')
BiocManager::install('pheatmap')
BiocManager::install('universalmotif')
BiocManager::install('seqLogo')
BiocManager::install('org.Mm.eg.db')
BiocManager::install('ATACseqQC')
BiocManager::install('JASPAR2020')
BiocManager::install('motifmatchr')
BiocManager::install('chromVAR')
BiocManager::install('ggseqlogo')
BiocManager::install('TFBSTools')
BiocManager::install('motifStack')
BiocManager::install('knitr')
BiocManager::install('testthat')
BiocManager::install('yaml')

# 查看是否都能加载
library('methods')
library('ggplot2')
library('rmarkdown')
library('ShortRead')
library('ashr')
library('ChIPQC')
library('DiffBind')
library('BSgenome.Hsapiens.UCSC.hg19')
library('Rsubread')
library('Rbowtie2')
library('R.utils')
library('Rsamtools')
library('BSgenome.Hsapiens.UCSC.hg38')
library('rtracklayer')
library('ChIPseeker')
library('soGGi')
library('GenomicAlignments')
library('TxDb.Hsapiens.UCSC.hg19.knownGene')
library('DESeq2')
library('BSgenome.Mmusculus.UCSC.mm10')
library('TxDb.Hsapiens.UCSC.hg38.knownGene')
library('tracktables')
library('clusterProfiler')
library('TxDb.Mmusculus.UCSC.mm10.knownGene')
library('devtools')
library('tidyr')
library('DT')
library('dplyr')
library('rGREAT')
library('MotifDb')
library('Biostrings')
library('GenomicRanges')
library('pheatmap')
library('universalmotif')
library('seqLogo')
library('org.Mm.eg.db')
library('ATACseqQC')
library('JASPAR2020')
library('motifmatchr')
library('chromVAR')
library('ggseqlogo')
library('TFBSTools')
library('motifStack')
library('knitr')
library('testthat')
library('yaml')

2、背景介绍

ATAC-seq使用转座酶在测序之前高效地切割可接近的DNA，提供了一种在全基因组范围内绘制可接近/开放染色质的方法。

与其他技术相比，ATAC-seq的优点有：

少量的组织样本即可（>10000个细胞）

实验周期短（~4hrs）

• DNaseseq：酶消化法从转录因子结合位点周围的开放染色质中提取信号

• MNaseseq：酶消化法提取代表核小体定位的信号

• ATACseq：使用转座酶，并提供了一种从单个样品的转录因子结合位点和核小体位置同时提取信号的方法

2.使用数据说明

使用的数据是技能书提供的，大家下载一下就行

find $CONDA_PREFIX -name "asperaweb_id_dsa.openssh"
/home/xiaocainiao/miniconda3/envs/ATAC/etc/asperaweb_id_dsa.openssh

key=/home/xiaocainiao/miniconda3/envs/ATAC/etc/asperaweb_id_dsa.openssh

ascp -v -QT -l 300m -P 33001 -k 1 -i $key era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR891/SRR891269/SRR891269_1.fastq.gz ./
ascp -v -QT -l 300m -P 33001 -k 1 -i $key era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR891/SRR891269/SRR891269_2.fastq.gz ./

3.参考序列数据

# human GRCh38
https://www.encodeproject.org/files/ENCFF356LFX/@@download/ENCFF356LFX.bed.gz
# 小鼠 mm10
https://www.encodeproject.org/files/ENCFF547MET/@@download/ENCFF547MET.bed.gz

4.ATACseq数据比对

4.1 创建参考基因组

## 4.1 创建参考基因组
library(BSgenome.Hsapiens.UCSC.hg19)
mainChromosomes <- paste0("chr",c(1:21,"X","Y","M"))
mainChrSeq <- lapply(mainChromosomes,
                     function(x)BSgenome.Hsapiens.UCSC.hg19[[x]])
names(mainChrSeq) <- mainChromosomes
mainChrSeqSet <- DNAStringSet(mainChrSeq)
writeXStringSet(mainChrSeqSet, "BSgenome.Hsapiens.UCSC.hg19.mainChrs.fa")

会在工作目录中生成一个fa文件：BSgenome.Hsapiens.UCSC.hg19.mainChrs.fa

4.2 创建 Rsubread 索引

## 4.2 创建 Rsubread 索引
library(Rsubread)
buildindex("BSgenome.Hsapiens.UCSC.hg19.mainChrs",
           "BSgenome.Hsapiens.UCSC.hg19.mainChrs.fa",
           indexSplit = TRUE,
           memory = 4000)

生成的index文件

4.3 比对 subread

# 比对
read1 <- "GSE47753/SRR891269_1.fastq.gz"
read2 <- "GSE47753/SRR891269_2.fastq.gz"

align("subread_index/BSgenome.Hsapiens.UCSC.hg19.mainChrs",
      readfile1=read1,
      readfile2=read2,
      output_file = "ATAC_50K_2.bam",
      nthreads=64,
      type=1,
      unique=TRUE,
      maxFragLength = 2000)

到这里就是生成比对结束的bam文件了，剩下还有作为sam、以及质控还有下游峰的一些分析