设备

此管道假定用户了解 shell 命令并且能够轻松地在基于 Linux 的操作系统上工作。

1.计算要求
以下过程可以在任何类似 Linux 的系统上运行。然而，该协议和公开可用的代码是为了通过由 SLURM 调度程序管理的高性能计算 (HPC) 集群执行命令而编写的。不过，此处提供的代码可以轻松转换为 TORQUE 或其他 HPC 系统。该管道假定工作 Perl 解释器版本 5.30.0 或更高版本以及 R 版本 3.6.2 或更高版本。

软件

以下分析过程使用了假定在用户 shell 环境中可用的几种标准计算工具：

1.BWA MEM （Li 和 Durbin，2009） v0.7.17； http://bio-bwa.sourceforge.net/bwa.shtml
2.SAM工具（李等人，2009） v1.14； http://www.htslib.org
3.床工具（昆兰和霍尔，2010） v2.27.1； https://bedtools.readthedocs.io/en/latest/
4.SRA 工具包（Leinonen 等人，2011） v2.11.1； https://github.com/ncbi/sra-tools
5.快速（陈等人，2018） v0.20.0； https://github.com/OpenGene/fastp
6.小猪v2.4； https://zlib.net/pigz/
7.MACS2 （刘，2014） v2.2.7.1； https://pypi.org/project/MACS2/
8.UCSC 二进制文件（Kent 等人，2010） v1.04.0； http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/
9.塔比克斯（李，2011） v0.2.6； http://www.htslib.org/doc/tabix.html
10.IGV （Thorvaldsdottir 等人，2013） v2.11.1； https://software.broadinstitute.org/software/igv
11.MEME （Grant 等人，2011） v5.4.1； https://meme-suite.org/meme/index.html
12.交叉映射（Zhao 等人，2014） v0.5.1； http://crossmap.sourceforge.net/
13.深度工具（拉米雷斯等人，2014） v3.5.1； https://deeptools.readthedocs.io/en/develop/index.html

输入数据

该计算管道的起始输入使用来自对玉米原生质体进行的 ATAC-STARR-seq 实验的双端测序数据（Ricci 等人， 2019）。 ATAC-STARR-seq 实验由 DNA 输入（ATAC-seq 文库）和 mRNA 读出（自转录调控区域）组成，用于识别表现出转录激活调控活性的基因组区域。

1.转染的 ATAC-seq DNA 输入 FASTQ
2.转录的 ATAC-seq mRNA FASTQ

程序

ATAC-STARR-seq 流程的概述如图2所示。

图2 .用于转座酶可及染色质分析以及自转录活性调控区测序 (ATAC-STARR-seq) 数据的计算工作流程。原始序列 ATAC-STARR-seq 数据首先从 NCBI GEO 获取，经过处理并使用BWA MEM与参考基因组比对，通过SAMtools过滤，使用BEDtools重新格式化为片段，并通过BEDtools进行比较和定制R 脚本为下游分析提供增强子活性的估计。

A.下载并准备必要的数据和参考基因组序列

1.玉米原生质体中转染玉米叶ATAC -seq 片段生成的原始 mRNA ATAC-STARR-seq 数据以及随附的 ATAC-seq 输入片段可在 NCBI GEO 上公开获得。 ATAC-STARR-seq mRNA 和 DNA 输入可以使用Fasterq -dump下载可从 SRA 工具包获取包裹：

#设置变量并下载 FASTQ 文件
目录 FASTQ_文件
cd FASTQ_files
fastq -dump -o B73_maize_DNA_ input.fastq SRR10964904
fastq -dump -o B73_maize_mRNA_output.fastq SRR10964905

2.pigz压缩 FASTQ 文件。默认情况下， pigz使用所有可用处理器，如果可用处理器数量未知，则使用 8 个处理器。或者，用户可以使用unix工具gzip来压缩 STARR mRNA 和 ATAC 输入 DNA FASTQ 文件。

# 压缩fastq文件
猪猪 *。快q

＃不跑
# 提示： gzip可以用作Pigz的替代品（并行gzip ）
#gzip _ *。快q

3.下载 B73 参考基因组序列和基因注释。原始文章将原始读数映射到 B73 玉米参考基因组的第 4 版。在本案例研究中，我将玉米 ATAC-STARR-seq 数据映射并分析到 B73 玉米参考基因组第 5 版，以展示更新的参考基因组和读取映射策略如何实现对公开数据集的信息性重新分析（Hufford 等人，2021））。

#下载参考数据
光盘.. /
目录基因组_参考
cd基因组_参考
获取 https://download.maizegdb.org/Zm-B73-REFERENCE-NAM-5.0/Zm-B73-REFERENCE-NAM-5.0.fa.gz
获取 https://download.maizegdb.org/Zm-B73-REFERENCE-NAM-5.0/Zm-B73-REFERENCE-NAM-5.0_Zm00001eb.1.gff3.gz

4.要将数据映射到参考基因组，我们首先需要解压缩 FASTA 文件。构建参考基因组索引是 BWA 比对的先决条件，并允许使用SAMtools和BEDtools等命令行工具更快地对 BAM/SAM/BED 格式文件进行比对后处理。

#为参考基因组 FASTA创建索引
gunzip Zm-B73-REFERENCE-NAM-5.0.fa.gz
萨姆工具 faidx Zm-B73-参考-NAM-5。 0.fa
bwa 索引 Zm-B73-REFERENCE-NAM-5。 0.fa

B.修剪适配器并删除低质量读数

1.如果 DNA 插入片段短于循环数，Illumina 平台可能会在 3' 端生成带有接头序列的读数。此外，由于定相，合成测序的保真度随着每个额外的循环而恶化，未去除的终止子帽导致的循环不同步最终导致后续循环中碱基识别的更大不确定性。去除接头污染和低质量碱基会增加可比对读数的总数，这在分析相对较低的序列复杂性实验（例如 ATAC-STARR-seq）时尤其重要。我们将使用fastp删除 mRNA 输出和 DNA 输入的测序适配器和低质量读数。可以在此处找到执行读取修剪的脚本： https: //github.com/Bio-protocol/Maize_ATAC_STARR_seq/blob/master/workflow/step01_trim_raw_reads.sh

#设置变量
光盘.. /
fastqdir =$PWD/ FASTQ_files
dna =B73_maize_DNA_input
rna =B73_maize_mRNA_output
线程=16

# 修剪和过滤 DNA 输入读数
fastp -j $ dna.json -h $dna.html -w $线程 \
-i $ fastqdir /${dna}_1.fastq.gz -I $ fastqdir /${dna}_2.fastq.gz \
-o $ fastqdir /${dna}_1.trim.fastq.gz -O $ fastqdir /${dna}_2.trim.fastq.gz

# 修剪和过滤 mRNA 输出读数
fastp -j $ rna.json -h $rna.html -w $线程 \
-i $ fastqdir /$ {rna} _1.fastq.gz -I $ fastqdir /$ {rna} _2.fastq.gz \
-o $ fastqdir /${rna} _1.trim.fastq.gz -O $ fastqdir /$ { rna } _2.trim.fastq.gzbS - | samtools排序 - > $ outdir /B73_maize_mRNA_output.raw.bam

#整理日志文件
目录 fastp_log_files
MV *. json *.html fastp_log_files

C.对齐和处理测序读取

1.经过修剪和质量过滤读数后，我们将输入 DNA 和输出 mRNA 读数与玉米 B73 版本 5 参考基因组进行比对。为了加速对齐和下游处理，我们使用 24 个 CPU (-t 24) 来对齐输入和输出读取。但是，用户应修改此值以反映系统上可用内核的数量。此外，我们将分割命中标记为二级比对（-M），以便在下游过滤掉，因为玉米基因组是高度重复的。 BWA MEM的输出通过管道传送到SAMtools 压缩视图（SAM 到 BAM）并在进一步处理之前按对齐坐标排序，以最大限度地减少磁盘占用空间。可以在此处找到执行这些步骤的脚本： https://github.com/Bio-protocol/Maize_ATAC_STARR_seq/blob/master/workflow/step02_align_STARR_data.sh

光盘.. /
目录 BAM_文件
outdir =$PWD/ BAM_files
refdir =$PWD/ Genome_Reference
fastqdir =$PWD/ FASTQ_files

#将DNA 输入和管道与samtools对齐，以进行 SAM 到 BAM 的转换
bwa mem -M -t 24 $ refdir /Zm-B73-REFERENCE-NAM-5。 0.fa $ fastqdir /B73_maize_DNA_input_1.trim.fastq.gz $ fastqdir /B73_maize_DNA_input_2.trim.fastq.gz | samtools视图 - bS - | samtools排序 - > $ outdir /B73_maize_DNA_input.raw.bam

#将RNA 输出和管道对齐到samtools以进行 SAM 到 BAM 的转换
bwa mem -M -t 24 $ refdir /Zm-B73-REFERENCE-NAM-5。 0.fa $ fastqdir /B73_maize_mRNA_output_1.trim.fastq.gz $ fastqdir /B73_maize_mRNA_output_2.trim.fastq.gz | samtools视图 - bS - | samtools排序 - > $ outdir /B73_maize_mRNA_output.raw.bam

2.为了确保只保留高质量的比对，这里我们删除不正确配对的读取（-f 3）、二次命中（-F 256）、低映射质量的比对（-q 10）和多个映射读取（ grep - v -E -e '\ bXA:Z :') 使用SAMtools和unix命令的组合。通过在SAMtools中设置 -h 标志，包含参考基因组和读取映射参数信息的标头保留在输出中查看命令。

# 过滤DNA输入比对
samtools view -h -q 10 -f 3 -F 256 $ outdir /B73_maize_DNA_input.raw.bam | grep -v -E -e '\ bXA:Z :' | samtools视图 - bSh -> $ outdir /B73_maize_DNA_input.mq10.pp.unique.bam

# 过滤 mRNA 输出比对
samtools view -h -q 10 -f 3 -F 256 $ outdir /B73_maize_ mRNA_output.raw.bam | grep -v -E -e '\ bXA:Z :' | samtools视图 - bSh -> $ outdir /B73_maize_mRNA_output.mq10.pp.unique.bam

3.ATAC-seq 和 STARR-seq 数据分析之间的主要区别在于如何通过测序捕获分析信息。对于配对末端 ATAC-seq，由于 Tn5 在其结合的基因组位置附近插入测序接头，因此染色质可及性被编码为测序读数的 5' 端。相比之下，STARR-seq 从能够激活自身转录的片段生成 mRNA 转录本；因此，整个 STARR-seq mRNA 和 DNA 片段可为分析提供信息。以下命令利用 BAM 双端比对中的 CIGAR 字符串来提取测序片段的坐标，以进行 DNA 输入和 mRNA 输出。可以在此处找到提取片段的脚本： https://github.com/Bio-protocol/Maize_ATAC_STARR_seq/blob/master/workflow/step03_extract_fragments.sh

#创建输出目录
目录 BED_文件

# 变量
outdir =$PWD/ BAM_files
beddir =$PWD/ BED_files
dna =B73_maize_DNA_input
rna =B73_maize_mRNA_output

bedtool的关于“缺失”配偶对的警告，这些警告反映了在先前步骤中过滤了配偶的一对）
echo“从 STARR DNA 输入中提取片段...”
samtools sort -n $ outdir /$dna.mq10.pp.unique.bam \
|床具 bamtobed - bedpe -我- \
|排序 -k1,1 -k2,2n - \
|剪切 -f1,2,6 -> $ beddir /$ dna.fragments.bed

# 提取mRNA输出片段
echo“从 STARR mRNA 输出中提取片段...”
samtools sort -n $ outdir /$rna.mq10.pp.unique.bam \
|床具 bamtobed - bedpe -我- \
|排序 -k1,1 -k2,2n - \
|切 -f1,2,6 -> $ beddir /$ rna.fragments.bed

D.识别背景下活动丰富的区域

1.为了识别具有激活转录能力的基因组区域，我们使用 MACS2 评估相对于输入 ATAC-seq 片段的 mRNA 读取的富集度。由于 MACS2 是通用峰值调用程序，因此我们需要调整默认设置以专门针对 ATAC-STARR-seq 数据定制分析。由于该实验没有使用独特的分子标识符，并且片段的转录本数量直接反映了其调节活性，因此保留了重复的 mRNA 片段（--keep-dup all）。为了直接使用由输入片段坐标确定的覆盖范围，我们关闭默认的片段移位模型（ --nomodel ）并将输入类型设置为BEDPE（-f BEDPE）。此外，我们通过将最小峰大小设置为 300 (--min-length 300) 来减小候选峰之间的最大间隙大小，以便在更广泛的监管区域 (--max-gap 50) 中识别精细映射的监管元素）。最后，我们使用背景覆盖率代替局部偏差，这有助于峰值检测（--nolambda ）。可以在此处找到执行峰值调用的脚本： https://github.com/Bio-protocol/Maize_ATAC_STARR_seq/blob/master/workflow/step04_call_peaks.sh

#准备输出目录和输入文件
目录峰值数据

#生成输入文件
uniq $ beddir /$ rna.fragments.bed > $ beddir /$ rna.fragments.uniq.bed
uniq $ beddir /$ dna.fragments.bed > $ beddir /$ dna.fragments.uniq.bed
猫 $ beddir /$ rna.fragments.uniq.bed $ beddir /$ dna.fragments.uniq.bed |排序-k1,1 -k2,2n -> $ beddir / ALL.fragments.uniq.bed

#查找监管区域
echo“调用监管区域而不重复删除......”
macs2 callpeak -t $ beddir /$ rna.fragments.bed \
-c $ beddir / ALL.fragments.uniq.bed \
--keep-dup 全部 \
--最大间隙 50 \
--最小长度 300 \
-- nolambda \
——无模型\
-f 贝德佩 \
-g 1.6e9 \
--bdg \ _
-n STARR_wdups

#查找监管区域
echo“调用重复删除的监管区域......”
macs2 callpeak -t $ beddir /$ rna.fragments.uniq.bed \
-c $ beddir /$ dna.fragments.uniq.bed \
--keep-dup 全部 \
--最大间隙 50 \
--最小长度 300 \
-- nolambda \
——无模型\
-f 贝德佩 \
-g 1.6e9 \
--bdg \ _
-n STARR_nodups

#使用所有独特片段查找调控区域
echo“通过聚合所有独特的片段来调用监管区域......”
macs2 callpeak -t $ beddir / ALL.fragments.uniq.bed \
--keep-dup 全部 \
--最大间隙 50 \
--最小长度 300 \
-- nolambda \
——无模型\
-f 贝德佩 \
-g 1.69e9 \
--bdg \ _
-n STARR_ALL

# 清理输出
mv STARR_*峰值数据

# 合并峰值
cd峰值数据
cat STARR_wdups_peaks.narrowPeak _ STARR_nodups_peaks.narrowPeak STARR_ALL_peaks.narrowPeak | STARR_ALL_peaks.narrowPeak |排序 -k1,1 -k2,2n - | bedtools合并 - i - > STARR_merged_peaks.bed

2.为了精细地估计监管活动，首先我们需要根据 mRNA 和 DNA 片段创建一个独特的区间列表。接下来，我们计算每个独特间隔的 mRNA 和 DNA 片段的交集。最后，我们删除所有临时文件以减少磁盘占用空间。可以在此处找到估计增强子活性的脚本： https: //github.com/Bio-protocol/Maize_ATAC_STARR_seq/blob/master/workflow/step05_estimate_enhancer_activity.sh

# 变量
beddir =$PWD/ BED_files
dna =B73_maize_DNA_input
rna =B73_maize_mRNA_output
ref=./ Genome_Reference /Zm-B73-REFERENCE-NAM-5.0.fa.fai

#对引用进行排序
排序 -k1,1 -k2,2n $ref > $ ref.sorted

# 合并RNA/DNA
猫 $ beddir /$ rna.fragments.uniq.bed $ beddir /$ dna.fragments.uniq.bed \
|排序 -k1,1 -k2,2n - \
|床具基因组cov - i - - bga -g $ ref.sorted \
|排序 -k1,1 -k2,2n - \
|切-f1-3-> $ beddir / Unique_genomic_intervals.bed

# 计算碎片数
bedtools intersect -a $ beddir / Unique_genomic_intervals.bed \
-b $ beddir /$ rna.fragments.bed \
-c -sorted -g $ ref.sorted > $ beddir /$ rna.activity.raw.bed

bedtools intersect -a $ beddir /$ rna.activity.raw.bed \
-b $ beddir /$ dna.fragments.bed \
-c -sorted -g $ ref.sorted > $ beddir /B73_maize_mRNA_DNA.activity.raw.bed

#清理临时文件
rm $ beddir / Unique_genomic_intervals.bed
rm $ beddir /$ rna.activity.raw.bed

3.我们和其他人将增强子活性定义为检测中 DNA 片段以独特片段间隔以 log 2 (mRNA/DNA) 形式产生的 mRNA 转录物的富集。在采用 mRNA 与 DNA 的 log 2比率之前，我们将输入和输出标准化为百万分之一，以考虑测序深度和复杂性的差异。将伪计数1添加到具有至少一个 RNA 或 DNA 片段的任何区间。也可以使用> Rscript从命令行运行以下代码使用以下脚本的Estimate_Enhancer_Activity.R ： https://github.com/Bio-protocol/Maize_ATAC_STARR_seq/blob/master/workflow/bin/Estimate_Enhancer_Activity.R

#打开一个交互式 R 会话来估计增强子活性
cd $床目录
右

# 加载数据
a <- read.table ("B73_maize_mRNA_DNA.activity.raw.bed")

#重新格式化
行名(a) <- 粘贴(a$V 1,a $V2,a$V3,sep="_")
a[ ,1:3] <- NULL
a <- as.matrix (a)
名称(a) <- c( "mRNA", "DNA")
a <- a[ rowSums (a )!= 0,]
一个 <- 一个 + 1

#标准化
a <- a %*% diag (x=1e6/ colSums (a))
名称(a) <- c( "mRNA", "DNA")
一个 <-作为。数据框（a）

#估计增强子活性
a$enhancer_activity <- log2( a$mRNA / a$DNA )

#重新格式化输出
rownames (a) <- gsub ( " scaf _", " scaf ", rownames (a))
df <- data.frame ( do.call ( rbind , strsplit ( rownames (a), "_")))
df$X1 <- gsub ( " scaf "," scaf _", as.character (df$X1))
mrna <- df
DNA <- df
df$X4 <- a$enhancer_activity
mrna$X4 <- a$mRNA
dna$X4 <- a$DNA

# 负面活动上限为 0
df$X4 <- ifelse ( df$X4 < 0, 0, df$X4)

#将增强剂活动 BEDGRAPH 文件保存到磁盘
write.table ( df , file="B73_maize.enhancer_activity.bdg",quote=F, row.names =F, col.names =F, sep ="\t")
write.table ( mrna , file="B73_maize.mRNA.bdg",quote=F, row.names =F, col.names =F, sep ="\t")
write.table ( dna , file="B73_maize.DNA.bdg",quote=F, row.names =F, col.names =F, sep ="\t")

#退出交互模式
q( )

4.为了可视化增强子活性以及任何给定位点的归一化 mRNA 和 DNA 片段，上一步中bedGraph格式的每百万覆盖率值可以转换为 bigwig 文件（使用 UCSC Utils 的bedGraphToBigWig ），以便使用集成基因组查看器 (IGV) 轻松可视化) 或JBrowse实例。

bedGraphToBigWig B73_maize.enhancer_activity.bdg ../ Genome_Reference /Zm-B73-REFERENCE-NAM-5.0.fa.fai.sorted B73_maize.enhancer_activity.bw

bedGraphToBigWig B73_maize.mRNA.bdg ../ Genome_Reference /Zm-B73-REFERENCE-NAM-5.0.fa.fai.sorted B73_maize.mRNA.bw

bedGraphToBigWig B73_maize.DNA.bdg ../ Genome_Reference /Zm-B73-REFERENCE-NAM-5.0.fa.fai.sorted B73_maize.DNA.bw

5.要查看增强子活性、mRNA 和 DNA 片段大佬文件，请在本地计算机上下载并安装 IGV ( https://software.broadinstitute.org/software/igv/download )。
6.解压、 bgzip并索引基因产物注释。然后，使用 IGV 中的File > Load from File…加载所有 bigwig、 narrowPeak和基因组注释文件。 IGV 屏幕截图的示例如图 3所示。

# 改变目录
cd ../基因组参考

# 解压
gunzip Zm-B73-REFERENCE-NAM-5.0_Zm00001eb.1.gff3.gz

# 按坐标排序并删除整个染色体间隔
排序-k1,1 -k4,4n Zm-B73-REFERENCE-NAM-5.0_Zm00001eb。 1.gff 3 | grep -v 程序集-> Zm-B73-REFERENCE-NAM-5.0_Zm00001eb.1.sorted.gff3

bgzip压缩
bgzip Zm-B73-REFERENCE-NAM-5.0_Zm00001eb。 1.排序.gff3

# 带有tabix 的索引
tabix -p gff Zm-B73-REFERENCE-NAM-5.0_Zm00001eb.1.sorted.gff3.gz

图3 .玉米中转座酶可及染色质分析分析的可视化，以及自转录活性调控区测序 (ATAC-STARR-seq) 数据。 DNA ATAC-seq 输入（蓝色）、自转录 mRNA 片段（粉色）和增强子活性的标准化（每百万读数）覆盖率 [紫色； log 2 (mRNA/DNA)] 23 kb 窗口。调节区域显示为黑条，而灰色环反映了预测的增强子基因链接。

E.创建控制区域列表

为了评估由 MACS2 确定的 STARR 调控区域相对于随机控制区域的富集度，我们首先需要识别在给定测序输出和读长的情况下可以唯一映射的基因组区域。尽管有许多方法可用于估计可映射性，但我将说明一种使用根据本实验的测序参数定制的合成读数的简单方法。首先，我们使用wgsim为 mRNA 和 DNA 输入生成相同数量、具有相同测序长度（36 个核苷酸）的随机读取对提供给SAMtools 的工具。然后，将合成读数重新映射并按照原始 STARR-seq 测序实验进行统一处理，以识别可唯一映射的区域。通过将随机控制区域选择限制为独特的可映射基因组间隔，我们确保下游比较不会因可映射性和重复组成而产生偏差。可以在此处找到构建控制区域的脚本： https: //github.com/Bio-protocol/Maize_ATAC_STARR_seq/blob/master/workflow/step06_create_control_regions.sh

#进入“ Genome_Reference ”目录
光盘。 /基因组_参考

# 估计读取次数
mRNA_counts =$( samtools view -c ./ BAM_files /B73_maize_mRNA_output.raw.bam)
DNA_counts =$( samtools view -c ./ BAM_files /B73_maize_DNA_input.raw.bam)

#使用SAMtools包中的wgsim 生成与 mRNA 输出匹配的模拟读取
wgsim -1 36 -2 36 -d 300 -N $ mRNA_counts Zm-B73-REFERENCE-NAM-5。 0.fa模拟_STARR_mRNA_r1.fq 模拟_STARR_mRNA_r2.fq

#生成与 DNA 输入匹配的模拟读取
wgsim -1 36 -2 36 -d 300 -N $ DNA_counts Zm-B73-REFERENCE-NAM-5。 0.fa模拟_STARR_DNA_r1.fq 模拟_STARR_DNA_r2.fq

#压缩合成的fastq文件
猪猪 *。弗克

1.使用与原始 STARR-seq 数据相同的管道重新映射合成读取。

# 变量
outdir =$PWD/ BAM_files
refdir =$PWD/ Genome_Reference
ref=$ refdir /Zm-B73-REFERENCE-NAM-5.0.fa.fai.sorted
fastqdir =$ refdir

#对齐合成 mRNA 输出并通过管道传输至samtools以进行 SAM 到 BAM 的转换
bwa mem -M -t 24 $ refdir /Zm-B73-REFERENCE-NAM-5。 0.fa \
$ fastqdir /simulated_STARR_mRNA_r1.fq.gz \
$ fastqdir /simualted_STARR_mRNA_r2.fq.gz \
| samtools视图 - bS - \
| samtools排序->$ outdir / simulated_STARR_mRNA.raw.bam

#将合成 DNA 输入和管道与samtools对齐，以进行 SAM 到 BAM 的转换
bwa mem -M -t 24 $ refdir /Zm-B73-REFERENCE-NAM-5。 0.fa \
$ fastqdir /simulated_STARR_DNA_r1.fq.gz \
$ fastqdir /simulated_STARR_DNA_r2.fq.gz \
| samtools视图 - bS - \
| samtools排序->$ outdir / simulated_STARR_DNA.raw.bam

2.使用原始 STARR-seq 管道处理和过滤读数。

# 过滤合成 mRNA 比对
echo“过滤合成的 STARR mRNA 比对...”
samtools view -h -q 10 -f 3 $ outdir / simulated_STARR_mRNA.raw.bam \
| grep -v -E -e '\ bXA:Z :' \
| samtools视图 - bSh -> $ outdir /simulated_STARR_mRNA.mq10.pp.unique.bam

# 过滤合成 DNA 比对
echo“过滤 STARR DNA 比对...”
samtools view -h -q 10 -f 3 $ outdir / simulated_STARR_DNA.raw.bam \
| grep -v -E -e '\ bXA:Z :' \
| samtools视图 - bSh -> $ outdir /simulated_STARR_DNA.mq10.pp.unique.bam

3.识别独特的可绘制地图的区域。

# 提取mRNA输出片段
echo“从模拟的 STARR mRNA 输出中提取片段......”
samtools sort -n $ outdir /simulated_STARR_mRNA.mq10.pp.unique.bam \
|床具 bamtobed - bedpe -我- \
|排序 -k1,1 -k2,2n - \
|切 -f1,2,6 -> $ refdir /simulated_STARR_mRNA.fragments.bed

# 提取 DNA 输入片段
echo“从模拟的 STARR DNA 输入中提取片段......”
samtools sort -n $ outdir /simulated_STARR_DNA.mq10.pp.unique.bam \
|床具 bamtobed - bedpe -我- \
|排序 -k1,1 -k2,2n - \
|切-f1,2,6 -> $ refdir /simulated_STARR_DNA.fragments.bed

#合并所有片段（这一步按坐标排序可能需要一段时间）
猫 $ refdir / simulated_STARR_mRNA.fragments.bed $ refdir / simulated_STARR_DNA.fragments.bed \
|排序 -k1,1 -k2,2n \
| bedtools合并-i ->$ refdir / mappable_genomic_regions.bed

4.仅使用可映射区域并排除 MACS2 输出的假定监管区域来构建控制区域。

#创建控件
峰值=$PWD/ Peak_data / STARR_merged_peaks.bed
床具洗牌 - i $peaks \
-g $ref \
-incl $ refdir / mappable_genomic_regions.bed \
-排除$峰值\
|排序-k1,1 -k2,2n -> $PWD/ Peak_data / STARR_CONTROL.bed

F.比较增强子活性

1.确定 MACS2 预测的增强子输出的增强子活性以及阴性对照区域。

#创建目录来包含分析
光盘.. /
mkdir 01_Peak_Analysis
cd 01_峰值_分析

# 将最大增强子活性映射到假定的调控区域（ wdups ）
bedtools地图 -a ../ Peak_data / STARR_merged_peaks.bed -b ../ BED_files /B73_maize.enhancer_activity.bdg -o max -c 4 > STARR_merged_peaks.enhancer_activity.bed

# 将最大增强子活性映射到控制
bedtools地图 -a ../ Peak_data / STARR_CONTROL.bed -b ../ BED_files /B73_maize.enhancer_activity.bdg -o max -c 4 > STARR_CONTROL.enhancer_activity.bed

2.为了去除具有与背景相似的增强子活性的调控区域，我们使用基于匹配的对照区域的经验错误发现率（ eFDR ）来过滤 STARR 调控区域。用户指定的eFDR阈值可识别控制集中去除 1-eFDR 百分比的控制区域的增强子活性值。在此示例中，我们将 FDR 设置为广泛使用的比率 0.05。以下代码执行并绘制eFDR过滤和增强器活动分布，并且可以使用> Rscript从命令行运行 eFDR_Filter_STARR_Peaks.R 。基于源自控制区域的eFDR阈值过滤 STARR 峰如图 4A所示。可以在此处找到以下代码的 R 脚本： https: //github.com/Bio-protocol/Maize_ATAC_STARR_seq/blob/master/workflow/bin/eFDR_Filter_STARR_Peaks.R

#启动交互式 R 会话
> R

# 加载库
图书馆（秤）

# 加载数据
starr <- read.table (" STARR_merged_peaks.enhancer_activity.bed ")
con <- read.table (” STARR_CONTROL.enhancer_activity.bed ”)

#将缺失值设置为0
星$V4[星$V4=='.'] <- 0
con$V4[con$V4=='.'] <- 0

#转换为数字
starr$V4 <- as.numeric ( as.character (starr$V4))
con$V4 <- as.numeric ( as.character (con$V4))

#获取经验阈值
罗斯福<- 0.05
阈值 <-分位数( con$V4, (1-fdr))

# 过滤 STARR 监管区域
过滤 <-子集( starr , starr$V4 >= 阈值)

#估计保留区域的比例
frac <-signif （ nrow （过滤）/ nrow （ starr ），数字=4）

#设置多面板绘图区域
pdf( “Density_eFDR_STARR_Peak_Filtering.pdf”，宽度=5，高度=5)

# 绘制具有重复项的 STARR 峰的对照/观察到的增强子活性
den.starr <- 密度(starr$V4)
den.con <- 密度(con$V4)
绘图（不适用，
xlab =“增强子活动”，
ylab =“密度”，
xlim =c(范围(范围( den.starr $x ), 范围( den.con $x))),
ylim =c(范围(范围( den.starr $y ), 范围( den.con$y ))))
多边形(x=c(min( den.rigid $x ), den.rigid$x , max( den.rigid$x )),
y= c( 0, den.starr$y , 0), col=alpha("darkorchid4", 0.5), border=NA)
多边形(x=c(min( den.con$x ), den.con$x , max( den.con$x )),
y= c( 0, den.con$y , 0), col=alpha("grey80", 0.5), border=NA)
abline ( v=阈值，col=“红色”， lty =2)
mtext ( paste0("STARR 峰通过 = ",frac," (", nrow (过滤), "/", nrow ( starr ),")"))

legend( "right", legend=c("STARR Peaks", "Control Peaks", Paste0(" eFDR = ", fdr )), col=c("darkorchid4", "grey75", "red"), border= c(NA, NA, “红色”), pch =c(16, 16, NA), lty =c(NA, NA, 2))

# 关闭设备
开发关闭( )

#保存过滤后的 STARR 监管区域
write.table （已过滤，文件=“STARR_merged_peaks.enhancer_activity.eFDR05.bed”，引用= F， row.names = F， col.names = F， sep =“\ t”）

# 退出R
q( )

图 4. 自转录活性调控区 (STARR) 调控区增强子活性分析。 A. STARR 峰（紫色）和随机对照区域（灰色）的增强子活性分布。红色虚线表示随机对照区域增强子活性的 95% 分位数。 B. ATAC-seq 输入（左）、mRNA 输出（中）和对照区域的增强子活性（右）、所有 STARR 峰和过滤的 STARR 的每百万读数 (RPM) 的平均（顶部）和单个位点热图峰值设置。

3.DeepTools评估过滤和未过滤 STARR 峰和对照的所有区域的相对增强子活性（图 4B ）。通过DeepTools绘制热图的脚本可以在这里找到： https://github.com/Bio-protocol/Maize_ATAC_STARR_seq/blob/master/workflow/step07_plot_enhancer_activity.sh

# 加载函数
获取地图( ){

＃输入
伊娜=1美元
id=$2
dat =.. /BED_files/*.bw

＃输出
outa=$id.mat.gz
outm=$id.mat.txt
图=$id.pdf

＃参数
线程=48
窗口=2000
仓=20
列= YlGnBu

#创建矩阵
计算矩阵参考点——参考点中心\
-S$数据\
-b $窗口 -a $窗口 \
-R $伊娜\
--缺失数据为零\
-o $outa \
-- outFileNameMatrix $ outm \
-p $threads -- binSize $bin

# 绘制热图
绘图热图——matrixFile $ outa \
-- colorMap YlGnBu \
-out $id.heatmap.pdf

}
导出-f getmaps

#为每个文件运行
获取地图 STARR_merged_peaks.enhancer_activity.bed STARR_峰值
getmaps STARR_merged_peaks.enhancer_activity.eFDR05.bed STARR_peaks_filtered
获取地图 STARR_CONTROL.enhancer_activity.bed 控制区域

G.玉米基因组中大调控域的鉴定

1.这些数据让我们提出一个问题：监管区域的大小与其增强子活动之间是否存在关系。哺乳动物系统中所谓的超级增强子描述了与细胞身份相关的高度活跃的转录激活调节域，与典型的增强子相比，其表现出增加的 TF 结合位点密度（赫尼兹等人。，2013）。 STARR 峰和增强子活性与其他数据集的整合使我们能够确定玉米中是否存在类似的过度活跃的调控域。为了查询 TF 结合位点密度，我们首先从MEME下载已知 TF 的位置权重矩阵数据库并使用fimo识别假定的 TF 结合位点（同样来自MEME suite)以小于 1 -5的P值阈值为条件。可以在此处找到识别大型监管域的脚本： https: //github.com/Bio-protocol/Maize_ATAC_STARR_seq/blob/master/workflow/step08_identify_large_regulatory_domains.sh

#为 TFBS 分析创建一个新目录
光盘.. /
mkdir 02_Hyperactive_Regulatory_Region_Analysis
光盘。 /02_Hyperactive_Regulatory_Region_Analysis

#下载并解压主题数据库
wget https://meme-suite.org/meme/meme-software/Databases/motifs/motif_databases.12.23.tgz
tar - xvzfmotive_databases.12.23.tgz
rm主题数据库.12.23.tgz

# 变量
线程=16
ref =.. / Genome_Reference /Zm-B73-REFERENCE-NAM-5.0.fa
峰=../01_Peak_Analysis/STARR_merged_peaks.enhancer_activity.eFDR05.bed
控制=../01_Peak_Analysis/ STARR_CONTROL.enhancer_activity.bed
主题=./主题数据库/ARABD/ArabidopsisDAPv1.meme

# 提取fasta序列
床具 getfasta -bed $peaks -fi $ref - fo $ peaks.fasta
床具 getfasta -bed $controls -fi $ref - fo $ controls.fasta

#识别假定的 TFBS
菲莫——OC _ TFBS_peaks $motifs $ peaks.fasta
菲莫——OC _ TFBS_controls $motifs $ controls.fasta

#使用github存储库中提供的perl脚本重新格式化fimo输出（过滤 p 值 > 1e-5）（ https://github.com/Bio-protocol/Maize_ATAC_STARR_seq/blob/master/workflow/bin/convertMotifCoord.pl）
perl ConvertMotifCoord.pl TFBS_peaks / fimo.gff | sed -e 's/_ tnt //g' - | sed -e 's/_ tnt //g' - |排序-k1,1 -k2,2n -> TFBS_peaks.motifs.bed
perl ConvertMotifCoord.pl TFBS_controls / fimo.gff | sed -e 's/_ tnt //g' - | sed -e 's/_ tnt //g' - |排序-k1,1 -k2,2n -> TFBS_controls.motifs.bed

#注释STARR 和对照峰的基序覆盖率/计数
bedtools注释 - i ../01_Peak_Analysis/STARR_merged_peaks.enhancer_activity.eFDR05.bed -files TFBS_peaks.motifs.bed -both |排序-k1,1 -k2,2n -> STARR_merged_peaks.enhancer_activity.eFDR05.ann.bed
bedtools注释 - i ../01_Peak_Analysis/ STARR_CONTROL.enhancer_activity.bed -files TFBS_controls.motifs.bed -both |排序-k1,1 -k2,2n -> STARR_CONTROL.enhancer_activity.ann.bed

# 提取基因
perl -ne ' if( $_ =~ /^#/){next;} chomp;my@col =split("\t",$_);if($col[2] eq '基因'){print “$_\n”;}' ../ Genome_Reference /Zm-B73-REFERENCE-NAM-5.0_Zm00001eb.1.gff3 |排序-k1,1 -k4,4n -> ../Genome_Reference/Zm-B73-REFERENCE-NAM-5.0_Zm00001eb.1.genes.gff3

#对 STARR 和对照峰的基因组背景进行分类（您可以忽略bedtools关于不一致命名约定的警告，您可以感谢基因组组装团队提供这些烦人但无害的警告）
bedtools最接近-a STARR_merged_peaks.enhancer_activity.eFDR05.ann.bed -b ../Genome_Reference/Zm-B73-REFERENCE-NAM-5.0_Zm00001eb.1.genes.gff3 -D b > STARR_merged_peaks.enhancer_activity.eFDR05.ann2.bed
bedtools最接近-a STARR_CONTROL.enhancer_activity.ann.bed -b ../Genome_Reference/Zm-B73-REFERENCE-NAM-5.0_Zm00001eb.1.genes.gff3 -D b > STARR_CONTROL.enhancer_activity.ann2.bed

＃清理_
mv STARR_merged_peaks.enhancer_activity.eFDR05.ann2.bed STARR_merged_peaks.enhancer_activity.eFDR05.ann.bed
mv STARR_CONTROL.enhancer_activity.ann2.bed STARR_CONTROL.enhancer_activity.ann.bed

2.我们现在可以研究调控区域大小、基序密度和增强子活性之间的关系，以确定假定的调控域（图 5A – 5F ）。为此，我们将启动交互式 R 会话并从上面加载带注释的峰值和控制文件。可以在此处找到自动执行以下代码的脚本： https: //github.com/Bio-protocol/Maize_ATAC_STARR_seq/blob/master/workflow/bin/Characterize_Regulatory_Regions.R

# 打开R
> R

# 分析监管区域

# 加载库
库（ viioplot ）
库（ dplyr ）
图书馆（大众）
库（ RColorBrewer ）
图书馆（秤）

# 加载数据
starr <- read.table ("STARR_merged_peaks.enhancer_activity.eFDR05.ann.bed")
控制 <- read.table (" STARR_CONTROL.enhancer_activity.ann.bed ")

#选择随机对照区域来匹配过滤后的 STARR 峰
控制 <- 控制[样本( nrow ( starr )),]

#为了清晰起见，重命名列（ frac_RR_motif = 基序覆盖的调控区域的分数）
斯塔尔[ ,7:15] <- NULL
控制[ ,7:15] <- NULL
colnames ( starr )[ 4:7] <- c("activity", " motive_counts ", " frac_RR_motif ", " gene_distance ")
colnames (对照)[ 4:7] <- c("activity", " motive_counts ", " frac_RR_motif ", " gene_distance ")

＃分类
starr$class <- ifelse ( ( starr$gene_distance < 0 & starr$gene_distance > -200), "TSS", ifelse ( starr$gene_distance < -200 & starr$gene_distance > -2000, "启动子", ifelse ( starr$基因距离> 0 & starr$gene_distance < 1000, "TTS", ifelse ( starr$gene_distance == 0, "基因", "基因间"))))

控制$class <- ifelse ( (控制$gene_distance < 0 & control$gene_distance > -200), "TSS", ifelse (控制$gene_distance < -200 &控制$gene_distance > -2000, "启动子", ifelse ( control$基因距离> 0 &控制$基因距离< 1000, "TTS", ifelse (控制$基因_距离== 0, "基因", "基因间"))))

# 绘图分布
pdf( “STARR_peak_control_genomic_distribution.pdf”，宽度=10，高度=5)
布局(矩阵( c( 1:2), nrow =1))
饼图（表（ starr$class ））
饼图（表（控制$类））
开发关闭( )

#估计监管区域大小（以 log10 为单位）
星$大小<- log10(星$V3-星$V2)
控制$大小<- log10(控制$V3-控制$V2)

#比较峰值和对照之间的大小（健全性检查）
pdf( “STARR_peak_control_sizes.pdf”，宽度=5，高度=6)
vioplot ( starr$size , control$size ,
ylab =“间隔大小（log10）”，
col= c( “ dodgerblue ”,“grey75”),
名称= c（ paste0（“STARR峰值\ n（n =”， nrow （ starr ），“）”），
Paste0("控制区域\n(n= ", nrow (控制),")")))
开发关闭( )

#比较主题数量
pval <-威尔科克斯。测试( starr$motif_counts , control$motif_counts )$ p.value
pval <- ifelse ( pval ==0, 2.2e-16, pval )
平均值.峰值<- 平均值( starr$motif_counts )
mean.cont <-mean( control$motif_counts )

#找到控制基序计数的 95% 分位数
upper.threshold <- 分位数( control$motif_counts , 0.95)

＃阴谋
pdf( “STARR_peak_control_motif_counts.pdf”，宽度=5，高度=6)
vioplot (log1p( starr$motif_counts ), log1p( control$motif_counts ),
ylab =“log2（基序计数 + 1）”，
col= c( “ dodgerblue ”,“grey75”),
名称= c（ paste0（“STARR峰值\ n（n =”， nrow （ starr ），“）”），
Paste0("控制区域\n(n= ", nrow (控制),")")),
ylim = c( 0,8),
面积等于T，
h=0.25）
mtext ( paste0("Wilcoxon Rank Sum P-value = ", signif ( pval ,digits=3)))
文本( 1, 7.5, labels=paste0("平均值 = ", signif (平均值.peak , 数字=3)))
文本( 2, 7.5, labels=paste0("平均值 = ", signif ( mean.cont , 数字=3)))
点( 1, log1p( upper.threshold ), col="red", pch ="-")
点( 2, log1p( upper.threshold ), col="red", pch ="-")
开发关闭( )

# 基于 95% 分位数控制距离按基序计数分割 STARR 区域
starr$group <- ifelse ( starr$motif_counts >= upper.threshold , "高", "低")
pval <-克鲁斯卡尔。测试( starr$activity , starr$group )$ p.value
pdf( “STARR_peak_activity_vs_group.pdf”，宽度=5，高度=6)
vioplot ( starr$activity~starr$group ，
ylab =“增强子活性”，
col= c( "dodgerblue4"," dodgerblue "),
名称= c（ paste0（“主题丰富\ n STARR峰值\ n（n =”， nrow （ starr [ starr $ group ==“high”，]），“）”），
Paste0("基序耗尽\n STARR 峰\n (n=", nrow ( starr [ starr$group =="low",]),")")),
面积等于=F,
xlab =“”，
h=0.25）
mtext ( paste0("Kruskal-Wallis 秩和 P 值 = ", signif ( pval ,digits=3)))
开发关闭( )

#比较STARR 区域大小
pval <-克鲁斯卡尔。测试( starr$size , starr$group )$ p.value
pval <- ifelse ( pval ==0, 2.2e-16, pval )
pdf( “STARR_peak_size_vs_group.pdf”，宽度=5，高度=6)
vioplot ( starr$size~starr$group ，
ylab =“间隔大小（log10）”，
col= c( "dodgerblue4"," dodgerblue "),
名称= c（ paste0（“主题丰富\ n STARR峰值\ n（n =”， nrow （ starr [ starr $ group ==“high”，]），“）”），
Paste0("基序耗尽\n STARR 峰\n (n=", nrow ( starr [ starr$group =="low",]),")")),
面积等于=F,
xlab =“”，
h=0.25）
mtext ( paste0("Kruskal-Wallis 秩和 P 值 = ", signif ( pval ,digits=3)))
开发关闭( )

#比较主题覆盖率
pval <-克鲁斯卡尔。测试( starr$frac_RR_motif , starr$group )$ p.value
pval <- ifelse ( pval ==0, 2.2e-16, pval )
pdf( “STARR_peak_motif_coverage_vs_group.pdf”，宽度=5，高度=6)
vioplot ( starr$frac_RR_motif~starr$group ，
ylab =“分数基序覆盖率”，
col= c( "dodgerblue4"," dodgerblue "),
名称= c（ paste0（“主题丰富\ n STARR峰值\ n（n =”， nrow （ starr [ starr $ group ==“high”，]），“）”），
Paste0("基序耗尽\n STARR 峰\n (n=", nrow ( starr [ starr$group =="low",]),")")),
面积等于=F,
xlab =“”）
mtext ( paste0("Kruskal-Wallis 秩和 P 值 = ", signif ( pval ,digits=3)))
开发关闭( )

#按组划分
starr.me <-子集( starr , starr$group =="high")
starr.md <-子集( starr , starr$group =="low")
write.table (starr.me, file="STARR_starrs_peaks.enhancer_activity.eFDR05.ann.high_motif.bed",
quote=F, row.names =F, col.names =F, sep ="\t")
write.table （starr.md，文件=“STARR_starrs_peaks.enhancer_activity.eFDR05.ann.low_motif.bed”，
quote=F, row.names =F, col.names =F, sep ="\t")

图5 .基序密集增强子调控域的鉴定。 A.自转录活性调节区 ( STARR) 峰（左）和对照区（右）的基因组分布。 B. 对照区域（灰色）和 STARR 峰（蓝色）间隔长度的分布。 C. 对照区域（灰色）和 STARR 峰（蓝色）中基序计数的分布。红色虚线表示来自控制区域的基序计数的 95% 分位数，用于将 STARR 峰分类为高基序计数类别和低基序计数类别。 D. 具有富集（深蓝色）和贫乏（浅蓝色）基序计数的 STARR 峰的增强子活性分布。 E. 基序富集（深蓝色）和基序耗尽（浅蓝色）STARR 峰的间隔长度分布。 F. 富含基序（深蓝色）和基序贫乏（浅蓝色）的 STARR 峰被基序覆盖的 STARR 峰分数的分布。 G. 热图显示了通过各种玉米细胞类型中相对染色质可及性缩放的基因间基序富集的 STARR 峰上的 Z 分数转换基序增强子活性。

3.为了确定大的基因间调控域区域是否与细胞身份相关，我们将利用来自多个玉米器官的最新单细胞 ATAC-seq ( scATAC -seq) 数据集，将增强子活性与各种细胞类型特异性 ACR 进行比较（马兰德等人， 2021）。首先，下载包含每种分析细胞类型的可访问染色质区域的标准化可访问性计数的矩阵。然后，我们提取 ACR 基因组坐标（位于 B73 参考基因组版本 4 中）并使用 CrossMap 将其转换为 B73 参考基因组版本5 工具和链文件。

#下载计数矩阵
wget -O maize_scATAC_atlas_ACR_celltype_CPM.txt.gz https://www.ncbi.nlm.nih.gov/geo/download/\?acc\=GSE155178\&format\=file\&file\=GSE155178%5Fmaize%5FscATAC%5Fatlas%5FACR %5Fcelltype%5FCPM%2Etxt%2Egz

# 解压
gunzip maize_scATAC_atlas_ACR_celltype_CPM.txt.gz

#下载链文件
wget https://download.maizegdb.org/Zm-B73-REFERENCE-NAM-5.0/chain_files/B73_RefGen_v4_to_Zm-B73-REFERENCE-NAM-5.0.chain

# 提取坐标并使染色体名称符合 V4 参考
切-f1 maize_scATAC_atlas_ACR_celltype_CPM.txt \
| grep '^chr' - \
| perl -ne ' chomp;my @col =split("_",$_); print"$col [0]\ t$col [1]\ t$col [2]\n";' -\
| sed -e 's/chrB73V4ctg/B73V4_ctg/g' - \
| sed -e 's/chr//g' - \
|排序-k1,1 -k2,2n -> maize_scATAC_atlas_ACRs.bed

#将ACR坐标从V4转换为V5
CrossMap.py 床 B73_RefGen_v4_to_Zm-B73-REFERENCE-NAM-5。 0.链 maize_scATAC_atlas_ACRs.bed > maize_scATAC_atlas_ACRs.V4_to_V5.txt

#丢弃未映射和分割的投影
grep -v '取消映射\|split' maize_scATAC_atlas_ACRs.V4_to_V5.txt \
| perl -ne 'chomp; my@col =split("\t ",$ _);打印“chr$col[0]_$col[1]_$col[2]\tchr$col[4]_$col[5]_$col[6]\n”;' -\
|排序 -k1,1 -k2,2n - \
| sed -e's/ chrscaf / scaf /g'-\
| sed -e 's/chrB73V4_ctg/chrB73V4ctg/g' - > maize_scATAC_atlas_ACRs.V4_to_V5.clean.txt

# 使用 R 更新矩阵文件中的 ACR 坐标
> R

#读入数据帧
conv <- read.table ("maize_scATAC_atlas_ACRs.V4_to_V5.clean.txt")
mat <- read.table ("maize_scATAC_atlas_ACR_celltype_CPM.txt")

# 投影后保留 ACR 的子集 mat 行
共享 <- intersect( rownames (mat), as.character (conv$V1))
垫 <- 垫[共享，]
行名称（转换）<- 转换 $V1
转换 <- 转换[共享,]

# 更新 mat rowID
行名(mat) <- conv$V2

#保存输出
write.table (mat, file="maize_scATAC_atlas_ACR_celltype_CPM.V5.txt", quote=F, row.names =T, col.names =T, sep = "\t")

# 退出R
q( )

#删除临时文件
rm maize_scATAC_atlas_ACR_celltype_CPM.txt maize_scATAC_atlas_ACRs.bed maize_scATAC_atlas_ACRs.V4_to_V5.txt maize_scATAC_atlas_ACRs.V4_to_V5.clean.txt

4.将scATAC -seq ACR 与具有富集基序计数的 STARR 峰相交。将scATAC -seq 矩阵和相交的 ACR/STARR 峰文件加载到 R 中，以估计按细胞类型的相对可及性缩放的增强子活性富集。由于 STARR-seq 数据源自玉米幼苗，我们进一步将scATAC -seq 细胞类型的分析限制于主要源自玉米幼苗的细胞类型。以下用 R 编写的代码可以使用名为“motif_enhancer_activity_maize_celltypes.R”的脚本执行，并提供按相对细胞类型可访问性缩放的各种基序的增强子活性估计，从而可以深入了解活性增强子的细胞类型特异性转录因子调节（图 5G ）。

# 提取ACR坐标
切-f1 maize_scATAC_atlas_ACR_celltype_CPM.V5.txt| grep -v '未知.5.50' | sed -e 's/ scaf _/ scaf /g' | sed -e 's/ scaf _/ scaf /g' | perl -ne ' chomp;my @col =split("_",$_); print"$col [0]\ t$col [1]\ t$col [2]\n";' - | sed -e 's/ scaf / scaf _/g' - |排序-k1,1 -k2,2n -> maize_scATAC_atlas_ACRs.V5.bed

#相交具有高基序数 STARR 峰的scATAC ACR
bedtools相交 -a STARR_starrs_peaks.enhancer_activity.eFDR05.ann.high_motif.bed -b maize_scATAC_atlas_ACRs.V5.bed - wa - wb > STARR_starrs_peaks.enhancer_activity.eFDR05.ann.high_motif.scATAC_ACRs.bed

# 将增强子活性映射到主题上
bedtools地图 -a TFBS_peaks.motifs.bed -b ../ BED_files /B73_maize.enhancer_activity.bdg -c 4 -o max > TFBS_peaks.motifs.enhancer_activity.bed

# 大型监管区域的主题
bedtools相交 -a TFBS_peaks.motifs.enhancer_activity.bed -b STARR_starrs_peaks.enhancer_activity.eFDR05.ann.high_motif.scATAC_ACRs.bed - wa - wb > TFBS_peaks.motifs.enhancer_activity.bed

# 打开 R （或者，可以在此处找到自动执行以下代码的脚本：https://github.com/Bio-protocol/Maize_ATAC_STARR_seq/blob/master/workflow/bin/motif_enhancer_activity_maize_celltypes.R）
> R

#估计增强子活性细胞类型特异性

# 加载库
库（ RColorBrewer ）
库（ gplots ）
库（ edgeR ）

# 加载数据
enh <- read.table ("STARR_starrs_peaks.enhancer_activity.eFDR05.ann.high_motif.scATAC_ACRs.bed")
acrs <- read.table ("maize_scATAC_atlas_ACR_celltype_CPM.V5.txt")
主题 <- read.table (" TFBS_peaks.motifs.ENRICHED.enhancer_activity.bed ")

# 代表性叶派生簇的子集
保持 <- c( "bulliform.2.26",
“束_鞘.2.16”，
“地面_meristem.7.69”，
“守卫_cell.7.74”，
“守卫_mother_cell.7.71”，
“L1_SAM.4.46”，
“叶_proangio.7.67”，
“叶肉.2.14”，
“薄壁组织.10.90”，
“原胚层.7.72”，
“气孔_前体.7.75”，
“子公司.7.68”）
all.acrs <- acrs

#重新缩放英亩
acrs <- cpm ( acrs , log=F)
acrs <- acrs [,保留]

# 按基因组特征划分的子集增强子
enh <-子集( enh , enh$V8=="基因间")

#从scATAC矩阵中获取重叠区域
enh$ids <- 粘贴(enh$V 11,enh $V12,enh$V13,sep="_")
enh <- enh [顺序（enh$V11，递减=T），]
enh <- enh [!重复( enh$ids ),]
共享 <- intersect( enh$ids , rownames ( acrs ))
行名( enh ) <- enh$ids
enh <- enh [共享，]
enh$starrIDs <-粘贴( enh$V1, enh$V2, enh$V3,sep="_")

# 过滤主题
图案$starrIDs <-粘贴（图案$V5，图案$V6，图案$V7， sep =“_”）
图案 <-图案[图案$starrIDs %in% 唯一( enh$starrIDs ),]

#标准化英亩
acrs <- t(应用( acrs , 1, 函数(x){x/max(x)}))

#迭代每个单元格类型
cts <-列名( acrs )
outs <- lapply ( cts , 函数(x){
访问 <- acrs [行名( enh ),x ]
名称（访问）<- enh$starrIDs
motifs.scores <- access[ motifs$starrIDs ] * as.numeric ( as.character (motifs$V15))
mtf <- data.frame (motif=motifs$V4, 分数= motive.scores )
aves <-聚合（分数~主题，数据= mtf ，FUN=平均值）
分数 <- aves$分数
名称(分数) <- aves$motif
返回（分数）
})
输出 <- do.call ( cbind , 输出)
colnames (出局) <- cts
vars <- apply( outs, 1, var)
输出 <-输出[变量 > 0,]
z <- as.matrix (t(scale(t(outs))))

# # 簇列
co <- hclust ( dist (t(outs) ))$顺序

# 重新排序行
z <- z [,co ]
row.o <- apply(z, 1, which.max )
z <- z[顺序( row.o , 递减=F),]

# 上限
z[ z < -3] <- -3
z[ z > 3] <- 3

#获取家人
tfs <- data.frame ( do.call ( rbind , strsplit ( rownames (z), "\\.")))
cols2 <- colorRampPalette ( brewer.pal (12, "配对"))(长度(唯一(tfs$X1)))
tfs$cols2 <- cols2[因子(tfs$X1)]

#可视化
pdf( “celltype_starr_motif_activity.pdf”，宽度=10，高度=10)
heatmap.2（z，scale =“无”，trace ='无'，
RowSideColors = tfs$cols ，
col= colorRampPalette ( rev( brewer.pal (9,“ RdBu ”)))(100),
useRaster =T， Colv =F， Rowv =F，树状图=“无”，边距= c（ 9,9））
开发关闭( )

致谢

这项研究得到了美国国家科学基金会 (DBI-1906869) 和国立卫生研究院 (1K99GM144742) 对 APM 的支持。本研究中分析的 ATAC-STARR-seq 数据最初由 Ricci、Lu、Ji 及其同事生成（里奇等人， 2019）。

利益争夺

APM 声明不存在竞争利益。

参考

Arnold, CD、Gerlach, D.、Stelzer, C.、Boryn, LM、Rath, M. 和 Stark, A. (2013)。 STARR-seq 鉴定的全基因组定量增强子活性图谱。科学339(6123)：1074-1077。
Buenrostro, JD、Giresi, PG、Zaba, LC、Chang, HY 和 Greenleaf, WJ (2013)。天然染色质转座，可快速、灵敏地对开放染色质、DNA 结合蛋白和核小体位置进行表观基因组分析。自然方法10(12)：1213-1218。
陈 S.、周 Y.、陈 Y. 和顾 J. (2018)。 fastp：一款超快的一体化 FASTQ 预处理器。生物信息学34(17)：i884-i890。
CE 格兰特、TL 贝利和 WS 诺布尔 (2011)。 FIMO：扫描给定主题的出现情况。生物信息学27(7): 1017-1018。
Hnisz, D.、Abraham, BJ、Lee, TI、Lau, A.、Saint-Andre, V.、Sigova, AA、Hoke, HA 和 Young, RA (2013)。控制细胞特性和疾病的超级增强剂。单元格155(4)：934-947。
Hufford, MB、Seetharam, AS、Woodhouse, MR、Chougule, KM、Ou, S.、Liu, J.、Ricci, WA、Guo, T.、Olson, A.、Qiu, Y. 等。（2021）。 26 个不同玉米基因组的从头组装、注释和比较分析。科学373(6555)：655-662。
Jores, T.、Tonnies, J.、Dorrity, MW、Cuperus, JT、Fields, S. 和 Queitsch, C. (2020)。通过 STARR-seq 鉴定烟叶中的植物增强剂及其组成元素[OPEN] 。植物细胞32(7)：2120-2131。
Kent, WJ、Zweig, AS、Barber, G.、Hinrichs, AS 和 Karolchik, D. (2010)。 BigWig 和 BigBed：支持浏览大型分布式数据集。生物信息学26(17)：2204-2207。
Leinonen, R.、Sukawara, H.、Shumway, M. 和国际核苷酸序列数据库 C. (2011)。顺序读取存档。核酸 Res 39（数据库问题）：D19-21。
李浩（2011）。 Tabix：从通用制表符分隔文件中快速检索序列特征。生物信息学27(5): 718-719。
Li, H. 和 Durbin, R. (2009)。通过 Burrows-Wheeler 变换快速准确地进行短读比对。生物信息学25(14)：1754-1760。
Li, H.、Handsaker, B.、Wysoker, A.、Fennell, T.、Ruan, J.、Homer, N.、Marth, G.、Abecasis, G.、Durbin, R. 和基因组计划数据处理， S.（2009）。序列比对/图谱格式和 SAM 工具。生物信息学25（16）：2078-2079。
刘，T.（2014）。使用基于模型的 ChIP-Seq 分析 (MACS) 来分析通过对胚胎干细胞中的蛋白质-DNA 相互作用进行测序而生成的短读段。方法分子生物学1150：81-95。
Marand, AP、Chen, Z.、Gallavotti, A. 和 Schmitz, RJ (2021)。单细胞分辨率的玉米顺式调控图谱。单元格184(11)：3041-3055。 e21。
Marand, AP、Zhang, T.、Zhu, B. 和 Jiang, J. (2017)。植物增强子的全基因组预测和表征。 Biochim Biophys Acta Gene Regul Mech 1860(1)：131-139。
Melnikov, A.、Murugan, A.、Zhang, X.、Tesileanu, T.、Wang, L.、Rogov, P.、Feizi, S.、Gnirke, A.、Callan, CG、Kinney, JB 等。（2012）。使用大规模并行报告分析对人类细胞中的诱导增强子进行系统剖析和优化。自然生物技术30(3)：271-277。
Minnoye, L.、Marinov, GK、Krausgruber, T.、Pan, L.、Marand, AP、Secchia, S.、Greenleaf, WJ、Furlong、EEM、Zhao, K.、Schmitz, RJ 等人。（2021）。染色质可及性分析方法。《自然评论方法入门》 1(1)：10。
AR 昆兰和 IM 霍尔 (2010)。 BEDTools：一套灵活的实用程序，用于比较基因组特征。生物信息学26(6): 841-842。
Ramirez, F.、Dundar, F.、Diehl, S.、Gruning, BA 和 Manke, T. (2014)。 deepTools：用于探索深度测序数据的灵活平台。核酸 Res 42（Web 服务器问题）：W187-191。
Ricci, WA、Lu, Z.、Ji, L.、Marand, AP、Ethridge, CL、Murphy, NG、Noshay, JM、Galli, M.、Mejia-Guerra, MK、Colome-Tatche, M. 等。（2019）。玉米基因组中广泛存在的长程顺式调控元件。自然植物5(12)：1237-1249。
Schmitz, RJ、Grotewold, E. 和 Stam, M. (2022)。植物中的顺式调控序列：其重要性、发现和未来的挑战。植物细胞34(2)：718-741。
孙静、何娜、牛丽、黄云、沉文、张云、李丽、侯成 (2019)。通过 STARR-seq 对水稻增强子进行全球定量作图。基因组学蛋白质组学生物信息学17(2): 140-153。
Thorvaldsdottir, H.、Robinson, JT 和 Mesirov, JP (2013)。 Integrative Genomics Viewer (IGV)：高性能基因组数据可视化和探索。简要生物信息14(2): 178-192。
赵 H.、孙 ZF、王 J.、黄 HJ、Kocher, JP 和王 LG (2014)。 CrossMap：基因组组装之间坐标转换的多功能工具。生物信息学30(7)：1006-1007。

补充资料

数据和代码可用性：所有数据和代码均已存入 GitHub： https: //github.com/Bio-protocol/Maize_ATAC_STARR_seq