软件升级虽然是一件值得高兴的是，但是代码变化太大却不是一件好消息。比如说Seurat，这个单细胞分析最常用的R包，它的2.x版本和3.x版本的变化就是翻天覆地。

为了能够重现别人的代码，你可能需要重装2.3.4版本的Seurat，官方提供了安装脚本

1

source("https://gist.githubusercontent.com/satijalab/beb9bb50dedc75ee023bd5d9be5fe684/raw/e103577735a2fba9da2ccca14ce1ac33e46c1bc4/old_seurat.R")

但是在window下安装会出一个报错, 提示无法在Windows下安装MacOS的版本

这个时候需要手动安装Seurat，注意，这里要求有Rtools才能进行编译

1

install.packages('Seurat', repos = 'https://satijalab.org/ran', type = "source")

可以设置install.packages里的参数lib.loc，让Seurat安装到其他的文件加，就不会替换原来的Seurat，同时用library加载的时候，也需要设置lib.loc.

上面是简单可行易操作的方法，唯一的问题是你不能同时加载两个Seurat版本，当然你也不会想这样子做，所以这是一个伪需求。

下面就是瞎折腾环节， 我要将Seurat的2.3.4版本单独搞出一个R包，Seurat2，这样子就可以同时加载这两个R包。

下载Seurat 2.3.4

1

wget https://satijalab.org/ran/src/contrib/Seurat_2.3.4.tar.gz

解压缩它

1
2

tar xf Seurat_2.3.4.tar.gz
cd Seurat

删除MD5文件，因为它会做文件检验

1

rm MD5

修改里面所有的Seurat替换成Seurat2, seurat替换成seurat2

1
2

find . -type f -print0 | xargs -0 sed -i "s/Seurat/Seurat2/g"
find . -type f -print0 | xargs -0 sed -i "s/seurat/seurat2/g"

这种无差别的替换会有一个问题，会把一些这类http://www.satijalab.org/seurat非代码信息中的seurat替换成seurat2，不过这并不影响实际函数的使用。

将R/seurat.R重名为R/seurat2.R

1

mv R/seurat.R R/seurat2.R

之后打包修改后的文件

1

tar -czf Seurat2.tar.gz Seurat

于是我们就可以安装我们自己修改后的R包了

1

install.packages("Seurat2.tar.gz", repos=NULL, type="source")

注意: 目前未修改测试数据集的对象，所以不能用Seurat2来运行pbmc_small的例子

由于我已经修改好了，所以你们也不需要自己再去修改了。我将其上传到GitHub，所以可以用devtools进行安装。

1

devtools::install_github("xuzhougeng/Seurat2")

由于代码只是简单修改，存在bug，其实不建议尝试使用。而且安装完Seurat2之后，我定性的评估了下Seurat2和Seurat3，发现两者在聚类分析上并没有太大区别, 所以更推荐大家使用最新的Seurat.

BigWig是一个能用于加载到基因组浏览器上展示的格式。它的格式比较复杂，不适合直接阅读，通常由BedGraph文件转换而来。

在R语言中可以通过rtracklayer的export.bw输入和输出BigWig文件。

默认情况下，我们导入的数据集是GRanges格式

1
2
3
4

test_path <- system.file("tests", package = "rtracklayer")
test_bw <- file.path(test_path, "test.bw")
gr <- import(test_bw)
gr 

输出信息如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

gr
GRanges object with 9 ranges and 1 metadata column:
      seqnames    ranges strand |     score
         <Rle> <IRanges>  <Rle> | <numeric>
  [1]     chr2     1-300      * |        -1
  [2]     chr2   301-600      * |     -0.75
  [3]     chr2   601-900      * |      -0.5
  [4]     chr2  901-1200      * |     -0.25
  [5]     chr2 1201-1500      * |         0
  [6]    chr19 1501-1800      * |      0.25
  [7]    chr19 1801-2100      * |       0.5
  [8]    chr19 2101-2400      * |      0.75
  [9]    chr19 2401-2700      * |         1
  -------
  seqinfo: 2 sequences from an unspecified genome

我们可以通过设定import里的as参数，使其读取为RleList

1
2

rle <- import(test_bw, as = "RleList")
rle

输出信息如下

1
2
3
4
5
6
7
8
9
10

RleList of length 2
$chr2
numeric-Rle of length 243199373 with 5 runs
  Lengths:       300       300       300       300 243198173
  Values :        -1     -0.75      -0.5     -0.25         0

$chr19
numeric-Rle of length 59128983 with 6 runs
  Lengths:     1500      300      300      300      300 59126283
  Values :        0     0.25      0.5     0.75        1        0

如果要输出BigWig文件，也只要保证提供的对象是GRanges或RleList， 只不过要保证一定要有score列表示信号强度。

1

export.bw(rle, "test.bw")

假如要用R语言输出下面这个GRanges对象为BigWig

1
2
3
4
5

gr <- GRanges(
        seqnames=Rle(c("chr1", "chr2", "chr1", "chr3"), c(1, 3, 2, 4)),
        ranges=IRanges(1:10, end=10),
        strand=Rle(strand(c("-", "+", "*", "+", "-")), c(1, 2, 2, 3, 2)),
        seqlengths=c(chr1=11, chr2=12, chr3=13))

那么思路就是

1. 获取染色体的长度
2. 将染色体按照一定宽度分窗, 可用函数为tile和slidingWindows
3. 统计每个窗口的read数, 需要用到函数coverage, Views和ViewSums
4. 输出BigWig

代码实现:

1
2
3
4
5
6
7
8
9
10
11
12
13
14

# 获取染色体长度
chromSizes <- GRanges(c("chr1","chr2","chr3"), IRanges(1,c(11,12,13)))
# 按照宽度2进行分窗
windows <- tile(chromSizes, width = 2)
# 统计每个窗口的read数
cvg <- coverage(gr)
cov_by_wnd <- Views(cvg, windows)
sum_by_wnd <- viewSums(cov_by_wnd)
for(i in seq_along(windows)){
  mcols(windows[[i]])['score'] <- sum_by_wnd[[i]]
}
# 输出
seqlengths(windows) <- c(11,12,13)
export.bw(unlist(windows), "tmp.bw")

最近(2019.8)在NBT上发表了一篇单细胞ATAC-seq的文章，题为 “Massively parallel single-cell chromatin landscapes of human immune cell development and intratumoral T cell exhaustion”。这篇文章由ATAC-seq的发明者所在的实验室用了10X Genomic 公司开发的sc-ATAC-seq技术分析而来。

文章一共测了20w个细胞，GEO编号是GSE129785, 里面是分析后的数据，而原始的FASTQ数据则在https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA532774, 数据量400G。此外，原始数据和处理过的数据还可以从https://support.10xgenomics.com/single-cell-atac/datasets进行下载。

我下载了其中66.7G的tar文件，解压后放在了百度网盘(链接：https://pan.baidu.com/s/1fu4MOhhW3aFtInsiJ5zBWg 提取码：6j7d ).

最为更重要的是，文章提供了核心分析部分用到的代码，代码地址为：https://github.com/GreenleafLab/10x-scATAC-2019

因此，我们就可以根据文章的数据和代码来学习10X sc-ATAC-seq的数据分析方法, 我们的目标是尽量重复出文章里出现的图，

文章数据

这篇文章能发到NBT，我觉得很大一部分的原因是它用大量的数据证明了scATAC-seq的可靠性，GSE129785编号中里放了67个样本，其中的GSE129785_RAW.tar(里面的文件是处理后的fragments.tsv.gz)有66.7G，解压后之后都是Cell Ranger ATAC (v 1.1.0) 处理后的fragments.gz文件, 我根据文章内容分门别类介绍这些数据：

第零类: 人类GM12878细胞系和小鼠A20细胞系(B lymphocytes)混合测序，通过滴定的方式来获取不同的细胞量，从而评估双/多细胞率(原文计算的结果是约1%)。一共有4个样本，如下

BY_500, BY_1k, BY_5k, BY_10k

第一类：人类PMBC中不同比例的不同类型细胞混合。作者的目的是为了分析多少细胞量会被鉴定出真实的类群，作者发现，百分之一或者千分之一的比例就可以鉴定出细胞类群，有14个样本

0p1_99p9_CD4Mem_CD8Naive, 0p1_99p9_Mono_T, 0p5_99p5_CD4Mem_CD8Naive, 0p5_99p5_Mono_T, 1_99_CD4Mem_CD8Naive, 1_99_Mono_T, 50_50_CD4Mem_CD8Naive, 50_50_Mono_T, 99_1_CD4Mem_CD8Naive, 99_1_MonoT, 99p5_0p5_CD4Mem_CD8Naive, 99p5_0p5_MonoT, 99p9_0p1_CD4Mem_CD8Naive, 99p9_0p1_Mono_T

其中，p是point的简写，0p1_99p9_CD4Mem_CD8Naive就是1 CD4Mem : 999 CD8Naive. 其中MonoT是作者手误，应该是Mono_T。文章里面只提到了monocytes 和 T cells，但是这里面的数据表明还有不同比例的CD细胞混合。

第二类：fresh PBMC和viably frozen PBMC，后者分为分选和未分选。作者的目的是为了确定即便是冻存的细胞也是能够做sc-ATAC-seq。有3个样本

Fresh_pbmc_5k, Frozen_sorted_pbmc_5k, Frozen_unsorted_pbmc_5k。

第三类：16个人健康人的peripheral blood, 4个重复

PBMC_Rep1, PBMC_Rep2, PBMC_Rep3, PBMC_Rep4

第四类：16个人健康人的 bone marrow cells, 就1个样本

Bone_Marrow_Rep1

第五类：16个人健康人的富集特定表面标记的细胞,

Dendritic_Cells, Monocytes, B_Cells, CD34_Progenitors_Rep1, Regulatory_T_Cells, Naive_CD4_T_Cells_Rep1, Memory_CD4_T_Cells_Rep1, CD4_HelperT, CD4_Memory, CD4_Naive, NK_Cells, Naive_CD8_T_Cells, Memory_CD8_T_Cells, Dendritic_all_cells, CD34_Progenitors_Rep2, Memory_CD4_T_Cells_Rep2, Naive_CD4_T_Cells_Rep2

第三，四，五类其实都是健康人的免疫细胞, 文章过滤后有61,806个细胞, 聚类结果见附录sup fig. 3a

第六类：免疫疗法前后的细胞，site-matched serial tumor biopsies pre- and post-PD-1 blockade (pembrolizumab) from five patients, plus post-therapy biopsies from two additional patients, 一共有21个样本，过滤后有37,818个细胞

• SU001_Immune_Post2, SU001_Tcell_Post2, SU001_Tcell_Post, SU001_Total_Post2, SU001_Total_Pre, SU001_Tumor_Immune_Post
• SU006_Immune_Pre, SU006_Tcell_Pre, SU006_Total_Post, SU006_Tumor_Pre
• SU008_Immune_Post, SU008_Immune_Pre, SU008_Tcell_Post, SU008_Tcell_Pre, SU008_Tumor_Post, SU008_Tumor_Pre
• SU009_Tcell_Post, SU009_Tcell_Pre, SU009_Tumor_Immune_Post, SU009_Tumor_Immune_Pre
• SU010_Total_Post, SU010_Total_Pre
• SU005_Total_Post
• SU007_Total_Post

注: 这里前5个病人还可以分类，一种是直接测Total，也就是unbiased fashion，另一种区分了Immune, Tcell, Tumor_Immune, 也就是cell sorting后分开测, T cells (CD45 + CD3 + ), non-T immune cells (CD45 + CD3 − ), stromal and tumor cells (CD45 − ), 参考附录sup fig. 6a

对于这些数据，文章用的是下面这套流程来进行处理

数据预处理

文章的数据预处理交给了Cell Ranger ATAC (v 1.1.0) ，用的是hg19 作为参考基因组，并不是最新GRCh38。

经@Jimmy @王诗翔 @伊现富老师 解释，是因为hg19的惯性太大，”周边”太多, 整个生信软件生态没有一致性更新，所以很多分析只能用hg19做。

Cell Ranger ATAC (v 1.1.0) 能做很多工作，

• Read过滤和比对(Read filtering and alignment)
• Barcode记数(Barcode counting)
• 鉴定转座酶切割位点(Identification of transposase cut sites)
• 发现可及染色质峰(Detection of accessible chromatin peaks)
• 定义细胞(Cell calling)
• peak和TF的count矩阵(Count matrix generation for peaks and transcription factors)
• 降维(Dimensionality reduction)
• 细胞聚类(Cell clustering)
• 聚类间差异可及性分析(Cluster differential accessibility)

原文用了很长一段描述(见Data processing using Cell Ranger ATAC software)来说明了Cell Ranger ATAC的具体过程。当然我们不在乎过程，我们只在乎结果。而结果就是下面这些:

• singlecell.csv: 每个细胞的信息，如是否在TSS
• possorted_bam.bam,possorted_bam.bam.bai: 处理后的BAM与其索引
• raw_peak_bc_matrix.h5: 原始的peak-cell矩阵, h5格式存放
• raw_peak_bc_matrix: 原始的peak-cell矩阵
• analysis: 各种分析结果，如聚类，富集，降维等
• filtered_peak_bc_matrix.h5: 过滤后的peak-cell矩阵, h5格式存放
• filtered_peak_bc_matrix: 过滤后的peak-cell矩阵
• fragments.tsv.gz, fragments.tsv.gz.tbi: 每个barcode的序列和它对应的基因组位置和数目
• filtered_tf_bc_matrix.h5: 过滤后的TF-cell矩阵, h5格式存放
• filtered_tf_bc_matrix: 过滤后的TF-cell矩阵
• cloupe.cloupe: Loupe Cell Browser的输入文件
• summary.csv, summary.json: 数据的统计结果, 以csv和json格式存放
• web_summary.html: 网页总结信息
• peaks.bed: 所有的peak汇总
• peak_annotation.csv: peak的注释结果

虽然结果很多，但是文章只用到了framgments.tsv.gz，就如同我们单细胞转录组分析也只要表达量矩阵。

scATAC-seq 数据分析

之后文章的方法部分就开始描述scATAC-seq的分析内容，一共有19项，描述多达3.5页。如果让作者在附录尽情发挥的话，那么每一项都会有1到2页内容，所以这将会一项非常持久的更新计划。

• Filtering cells by TSS enrichment and unique fragments
• Generating a counts matrix
• Generating union peak sets with LSI
• Reads-in-peaks-normalized bigwigs and sequencing tracks
• ATAC-seq-centric LSI clustering and visualization
• Inferring copy number amplification
• TF footprinting
• ChromVAR
• Computing gene activity scores using Cicero co-accessibility
• Analysis of autoimmune variants using Cicero co-accessibility and chromVAR.
• HiChIP meta-virtual 4C (metav4C) analysis for Cicero co-accessibility links
• Overlap of Cicero co-accessibility links with GTEx eQTLs
• Constructing ATAC-seq pseudo-bulk replicates of maximal variance
• Constructing gene score pseudo-bulk replicates of maximal variance
• Identification of cluster-specific peaks and gene scores through feature binarization
• Pseudotime analysis
• Barnyard mixing analysis
• Analysis of fresh versus frozen PBMCs
• Spike-in analysis

根据TSS富集和唯一片段过滤细胞

这一步用到的脚本是01_Filter_Cells_v2.R，读取的fragments为四列。前三列是fragment在染色体的位置，第四列是barcode信息，第五列则是该fragments被测了几次。

1
2

#chr	start	end	barcode	number
chr1    10059   10114   TTATGTCAGGTTGTTC-1      7

分析者先过滤不到那些fragments少于100的细胞。

之后，按照文章里提到过滤低质量细胞的标准进行过滤

• cut-offs of 1,000 unique nuclear fragments per cell
• a transcription start site (TSS) enrichment score of 8

每个细胞的唯一片段数计算比较简单也很好理解，而TSS enrichment score需要看他们之前在2018年发表在science文章的附录(Corces et. al 2018 science)

To get the fragment length distribution, the width of each
fragment/GRange was plotted. To get the TSS enrichment profile, each TSS from the R package “TxDb.Hsapiens.UCSC.hg38.knownGene” (accessed by transcripts(TxDb)) was extended 2000 bp in each direction and overlapped with the insertions (each end of a fragment) using “findOverlaps”. Next, the distance between the insertions and the strand-corrected TSS was calculated and the number of insertions occurring in each single-base bin was summed. To normalize this value to the local background, the accessibility at each position +/- 2000 bp from the TSS was normalized to the mean of the accessibility at positions +/-1900-2000 bp from the TSS. The final TSS enrichment reported was the maximum enrichment value within +/- 50 bp of the TSS after smoothing with a rolling mean every 51 bp.

这依赖于一个核心函数，insertionProfileSingles, 作用是分析每个细胞的插入谱

最后分析者为了避免潜在的doublets，还过滤了含有大于45,000的唯一fragmeng的细胞

构建counts矩阵

这部分描述包括后续的各种矩阵产生，包括cell-by-window, cell-by-peak, cell-by-tf。主要涉及02_Get_Peak_Set_hg19_v2.R 里的countInsertions函数。

第一步:获取 fragment GenomicRanges(数据结构如下)，也就是01_Filter_Cells_v2.R过滤后的结果

1
2
3
4
5
6
7
8
9

GRanges object with 6 ranges and 2 metadata columns:
      seqnames      ranges strand |                 RG         N
         <Rle>   <IRanges>  <Rle> |        <character> <numeric>
  [1]     chr1 10060-10114      * | TTATGTCAGGTTGTTC-1         7
  [2]     chr1 10060-10126      * | GAGGTCCCAGCGTCGT-1         4
  [3]     chr1 10061-10126      * | AACGTACAGTCTCTAG-1         5
  [4]     chr1 10061-10529      * | TAGCCGGAGGATGTAT-1         1
  [5]     chr1 10062-10114      * | TGAATCGTCTGACTTT-1        14
  [6]     chr1 10062-10126      * | CCGTAGGGTCACAGTT-1        13

第二步: 然后通过将fragments的起始位置和结束位置进行合并, 就得到了Tn5 insertion GenomicRange

1
2
3
4
5

by <- "RG"
inserts <- c(
    GRanges(seqnames = seqnames(fragments), ranges = IRanges(start(fragments), start(fragments)), RG = mcols(fragments)[,by]),
    GRanges(seqnames = seqnames(fragments), ranges = IRanges(end(fragments), end(fragments)), RG = mcols(fragments)[,by])
  )

第三步: 用findOverlaps找和insertion重叠的feature(如window, TSS, peak, TF)

1
2
3
4
5

overlapDF <- DataFrame(findOverlaps(query, 
                                    inserts, 
                                    ignore.strand = TRUE, 
                                    maxgap=-1L, minoverlap=0L, type = "any"))

第四步: 增加新的一列(id)到重叠输出对象中，这一列的内容是细胞barcode的ID，然后转成sparseMatrix提高存储和计算效率

1
2
3
4
5
6
7
8
9
10

overlapDF$name <- mcols(inserts)[overlapDF[, 2], by]
# transform类似于dplyr:mutate
overlapTDF <- transform(overlapDF, id = match(name, unique(name)))
# 构建稀疏矩阵对象
sparseM <- Matrix::sparseMatrix(
    i = overlapTDF[, 1], 
    j = overlapTDF[, 4],
    x = rep(1, nrow(overlapTDF)), 
    dims = c(length(query), length(unique(overlapDF$name))))
colnames(sparseM) <- unique(overlapDF$name)

第五步: Tn5插入数在peak的比例 (frip)计算方法为: 稀疏矩阵的列和和每个细胞的fragment数相除

1
2
3
4

inPeaks <- table(overlapDF$name)
total <- table(mcols(inserts)[, by])
total <- total[names(inPeaks)]
frip <- inPeaks / total

第六步(可选): counts矩阵用edgeR的cpm(matrix, log = TRUE, prior.count = 3)进行log标准化。这个矩阵用于降低那些低count值元素对变异的影响。该标准化假设不同细胞类型中所有染色质开放区的差异是比较小的。

1

logMat <- edgeR::cpm(clusterSums, log = TRUE, prior.count = 3)

利用LSI生成peak合集

方法参考自Cusanovich, D. A. et al (Cell 2018), 代码见02_Get_Peak_Set_hg19_v2.R, 步骤如下

第一步: 利用tile函数，沿着基因组构建2.5kb的window

1
2
3
4

genome <- BSgenome.Hsapiens.UCSC.hg19
chromSizes <- GRanges(names(seqlengths(genome)), IRanges(1, seqlengths(genome)))
chromSizes <- GenomeInfoDb::keepStandardChromosomes(chromSizes, pruning.mode = "coarse")
windows <- unlist(tile(chromSizes, width = 2500))

第二步: 利用countInsertions基于windows和fragments获取cell-by-window矩阵。

1
2
3
4
5
6
7
8
9

# lapply遍历所有的fragments, 然后合并
fragmentFiles <- list.files("data", pattern = ".rds", full.names = TRUE)
countsList <- lapply(seq_along(fragmentFiles), function(i){
	message(sprintf("%s of %s", i, length(fragmentFiles)))
	counts <- countInsertions(windows, readRDS(fragmentFiles[i]), by = "RG")[[1]]
  counts
})
mat <- lapply(countsList, function(x) x) %>% Reduce("cbind",.)
remove(countsList)

第三步: 矩阵二项化，并且只要所有细胞中前20,000最开放位点

1
2
3

# seuratLSI
mat@x[mat@x > 0] <- 1 
mat <- mat[head(order(Matrix::rowSums(mat),decreasing = TRUE),nFeatures),] 

第四步: 利用TF-IDF(term frequency-inverse document frequency)降维

1
2
3
4
5
6

# 计算单个细胞中每个开放位点的频率
freqs <- t(t(mat)/Matrix::colSums(mat))
# 计算逆文档频率
idf   <- as(log(1 + ncol(mat) / Matrix::rowSums(mat)), "sparseVector")
# 计算TF-IDF = TF * IDF
tfidf <- as(Matrix::Diagonal(x=as.vector(idf)), "sparseMatrix") %*% freqs

TF-IDF的思想: 如果某个词比较少见，但是它在这篇文章中多次出现，那么它很可能就反映了这篇文章的特性，正是我们所需要的关键词 – http://www.ruanyifeng.com/blog/2013/03/tf-idf.html

第五步: 利用R语言中的irlba对标准化TF-IDF进行奇异值分解(SVD) , 保留第2-25个维度(第一个维度通常和cell read depth高度相关)作为Seurat的输入

1
2
3
4
5
6
7
8
9
10
11

# SVD降维
svd <- irlba::irlba(tfidf, nComponents, nComponents)
svdDiag <- matrix(0, nrow=nComponents, ncol=nComponents)
diag(svdDiag) <- svd$d
matSVD <- t(svdDiag %*% t(svd$v))
rownames(matSVD) <- colnames(mat)
colnames(matSVD) <- paste0("PC",seq_len(ncol(matSVD)))
# 构建Seurat对象
obj <- CreateSeuratObject(mat, project='scATAC', min.cells=0, min.genes=0)
obj <- SetDimReduction(object = obj, reduction.type = "pca", slot = "cell.embeddings", new.data = matSVD)
obj <- SetDimReduction(object = obj, reduction.type = "pca", slot = "key", new.data = "PC")

第三，四，五步被整合到一个函数seuratLSI中， 文章用的是Seurat V2.3

第六步: 并用FindClusters进行SNN图聚类(默认0.8分辨率), 如果最小的细胞类群细胞数不够200，降低分辨率重新聚类， 一个函数addClusters实现。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

addClusters <- function(obj, minGroupSize = 50, dims.use = seq_len(50), initialResolution = 0.8){
  # 初步聚类
  currentResolution <- initialResolution
  obj <- FindClusters(object = obj, reduction.type = "pca", dims.use = dims.use, resolution = currentResolution, print.output = FALSE)
  # 获取最少聚类的细胞数和聚类数
  minSize <- min(table(obj@meta.data[[paste0("res.",currentResolution)]]))
  nClust <- length(unique(paste0(obj@meta.data[[paste0("res.",currentResolution)]])))
  message(sprintf("Current Resolution = %s, No of Clusters = %s, Minimum Cluster Size = %s", currentResolution, nClust, minSize))
  # 判断最少聚类细胞数是否大于200，如果低于200，降低分辨率(0.8 × N)
  while(minSize <= minGroupSize){
    obj@meta.data <- obj@meta.data[,-which(colnames(obj@meta.data)==paste0("res.",currentResolution))]
    currentResolution <- currentResolution*initialResolution
    obj <- FindClusters(object = obj, reduction.type = "pca", dims.use = dims.use, resolution = currentResolution, print.output = FALSE, force.recalc = TRUE)
    minSize <- min(table(obj@meta.data[[paste0("res.",currentResolution)]]))
    nClust <- length(unique(paste0(obj@meta.data[[paste0("res.",currentResolution)]])))
    message(sprintf("Current Resolution = %s, No of Clusters = %s, Minimum Cluster Size = %s", currentResolution, nClust, minSize))
  }
  return(obj)
}

需要注意，文章用的Seurat版本V2.3, 所以如果安装3.0版本, 上面的代码就失效了但是有解决方案，参考如下博客

• https://divingintogeneticsandgenomics.rbind.io/post/clustering-scatacseq-data-the-tf-idf-way/

第七步: 之后用MACS2对每个聚类里 Tn5-corrected singlebase insertions (each end of the Tn5-corrected fragments) 都进行peak calling，参数设置为–shift -75–extsize 150–nomodel–callsummits–nolambda–keep-dup all -q 0.05

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

# 保存BED文件
dirClusters <- "results/LSI-Cluster-Beds/"
dir.create(dirClusters)
for(i in seq_along(fragmentFiles)){
	fragments <-readRDS(fragmentFiles[i])
	for(j in seq_along(clusterResults)){
	  message(sprintf("%s of %s", j, length(clusterResults)))
	  fragmentsj <- fragments[fragments$RG %in% clusterResults[[j]]]
	  if(length(fragmentsj) > 0){
	    out <- data.frame(
	      chr = c(seqnames(fragmentsj), seqnames(fragmentsj)), 
	      start = c(as.integer(start(fragmentsj) - 1), as.integer(end(fragmentsj) - 1)), 
	      end = c(as.integer(start(fragmentsj)), as.integer(end(fragmentsj)))
	      ) %>% readr::write_tsv(
          x = ., 
          append = TRUE, 
          path = paste0(dirClusters, paste0(names(clusterResults)[j], ".bed")), 
          col_names = FALSE)
	    }
	}
}
# 运行MACS2
dirPeaks <- "results/LSI-Cluster-Peaks/"
method <- "q"
cutoff <- 0.05
shift <- -75
extsize <- 150
genome_size <- 2.7e9
for(j in seq_along(clusterResults)){
	message(sprintf("%s of %s", j, length(clusterResults)))
	clusterBedj <- paste0(dirClusters,names(clusterResults)[j],".bed")
	cmdPeaks <- sprintf(
	    "macs2 callpeak -g %s --name %s --treatment %s --outdir %s --format BED --nomodel --call-summits --nolambda --keep-dup all", 
	    genome_size, 
	    names(clusterResults)[j], 
	    clusterBedj, 
	    dirPeaks
	  )
	if (!is.null(shift) & !is.null(extsize)) {
	  cmdPeaks <- sprintf("%s --shift %s --extsize %s", cmdPeaks, shift, extsize)
	}
	if (tolower(method) == "p") {
	  cmdPeaks <- sprintf("%s -p %s", cmdPeaks, cutoff)
	}else {
	  cmdPeaks <- sprintf("%s -q %s", cmdPeaks, cutoff)
	}
	message("Running Macs2...")
	message(cmdPeaks)
	system(cmdPeaks, intern = TRUE)
}

之后读取peak summits文件，前后延伸250 bp, 最终宽度为501bp， 并基于ENCODE hg19 blacklist(https://www.encodeproject.org/annotations/ENCSR636HFF/)过滤，以及那些延伸后会超出染色体的peak。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

# 读取
df <- data.frame(
  samples = gsub("\\_summits.bed","",list.files(dirPeaks, pattern = "\\_summits.bed", full.names = FALSE)),
  groups = "scATAC",
  summits = list.files(dirPeaks, pattern = "\\_summits.bed", full.names = TRUE)
  )
# 合并peak, 延伸, 过滤
unionPeaks <- extendedPeakSet(
    df = df,
    BSgenome = genome, 
    extend = 250,
    blacklist = "data/hg19.blacklist.bed",
    nSummits = 200000
  )
# 过滤Y染色体和线粒体
unionPeaks <- unionPeaks[seqnames(unionPeaks) %in% paste0("chr",c(1:22,"X"))]
unionPeaks <- keepSeqlevels(unionPeaks, paste0("chr",c(1:22,"X")))

其中extendPeakSet里的函数nonOverlappingGRanges 会对单个样本中相互重叠的peak进行合并处理。步骤如下

• 如果存在有直接重叠的peak，就只保留最显著的peak，过滤不怎么显著的peak
• 迭代上一步，处理下一个重叠的peak，直到这类直接性重叠的peak要么都保留，要么都移除。这一步用于处理每个聚类的peak集，保留前200,000延伸的summit, 得到了每个cluster的’cluster-specific’ peak集
• 接着对每个样本的MACS2 peak score进行标准化，即−log10(Q value), 然后通过将每个个体的score转成中位数trunc(rank(v))/length(v)’, 其中v表示MACS2 peak scores 向量, 等到”score quantitle”。标准化之后peak就能够在聚类之间进行相互比较，能够产生每个数据集的peak合集
• 合并所有的cluster的peak，使用之前的迭代方法处理重叠，得到了最显著的peak，抛弃了不怎么显著的peak
• 接着将那些和基因组gap区(也就是序列是N)重叠peak移掉，以及Y染色体上的peak

第八步: 构建cell-by-peak矩阵. 用到的是countInsertions函数，只不过第一个参数改为了unionPeaks.

1
2
3
4
5
6

#Create Counts list
countsPeaksList <- lapply(seq_along(fragmentFiles), function(i){
  message(sprintf("%s of %s", i, length(fragmentFiles)))
  gc()
  countInsertions(unionPeaks, readRDS(fragmentFiles[i]), by = "RG")
})

做一些总结性分析

1
2
3
4
5
6
7
8
9
10
11
12

#CountsMatrix
mat <- lapply(countsPeaksList, function(x) x[[1]]) %>% Reduce("cbind",.)
frip <- lapply(countsPeaksList, function(x) x[[2]]) %>% unlist
total <- lapply(countsPeaksList, function(x) x[[3]]) %>% unlist

se <- SummarizedExperiment(
  assays = SimpleList(counts = mat), 
  rowRanges = unionPeaks
  )
rownames(se) <- paste(seqnames(se),start(se),end(se),sep="_")
colData(se)$FRIP <- frip
colData(se)$uniqueFrags <- total / 2

Reads-in-peaks-normalized bigwigs and sequencing tracks.

这部分内容在GitHub上没有代码，所以我根据文章描述编写代码实现

…

ATAC-seq-centric LSI 聚类和可视化

以Cusanovich et. al (Nature 2018) 提出的策略计算TF-IDF转换, 然后irlba进行奇异值分解(SVD), 接着用前50维作为Seurat的输入，实现SNN图聚类，然后用FindCluster寻找聚类。 你会发现这个内容和上一节类似，只不过上一节是基于window, 这一节是基于peak。

1
2
3
4
5
6

# TF-IDF, SVD 
obj <- seuratLSI(assay(se), nComponents = max(nPCs1), nFeatures = NULL)
obj@meta.data <- as.data.frame(cbind(obj@meta.data, colData(se)))
# 聚类分析
obj <- FindClusters(object = obj, dims.use = nPCs1, print.output = TRUE, n.start = 10)

分析者发现这里存在能被检测到的批次效应，它会混淆后续的分析。因此为了减缓影响，他们从二项化后到可及矩阵中计算每个聚类的和，然后对edgeR的cpm(matrix, log = TRUE, prior.count = 3)结果进行log标准化

1
2
3
4

mat <- assay(se)
mat@x[mat@x > 0] <- 1
clusterSums <- groupSums(mat = mat, groups = paste0("C",obj@meta.data$res.0.8), sparse = TRUE)
logMat <- edgeR::cpm(clusterSums, log = TRUE, prior.count = 3)

然后用rowVars挑选前25,000个波动peak，这是基于上一步标准化的矩阵

1

varPeaks <- head(order(matrixStats::rowVars(logMat), decreasing = TRUE), nTop)

为什么不用原来的稀疏二项矩阵，而用log标准化的CPM矩阵，原因有亮点

• 能够降低不同聚类间因细胞数不同导致的误差
• count转换到标准化空间后的均值-变异关系会减弱

之后重新用25,000波动peak提取稀疏的二项化矩阵，然后重新做TF-IDF, SVD, SNN graph聚类等运算。

1
2
3
4
5
6
7
8
9

#Re-run Seurat LSI
message("Making Seurat LSI Object...")
obj2 <- seuratLSI(assay(se)[varPeaks,], nComponents = max(nPCs2), nFeatures = NULL)
stopifnot(identical(rownames(obj2@meta.data), colnames(se)))
obj2@meta.data <- as.data.frame(cbind(obj2@meta.data, colData(se)))

message("Adding Graph Clusters...")
obj2 <- FindClusters(object = obj2, reduction.type = "pca", dims.use = nPCs2, print.output = TRUE, n.start = 10)

用Seurat::RunUMAP函数进行可视化展示

1
2
3
4
5
6
7
8
9
10
11
12
13

message("Running UMAP")
obj2 <- RunUMAP(object = obj2, reduction.use = "pca", dims.use = nPCs2)
plotUMAP <- data.frame(GetCellEmbeddings(obj2,reduction.type="umap"), obj2@meta.data)
colnames(plotUMAP) <- c("x","y",colnames(plotUMAP)[3:ncol(plotUMAP)])
clustCol <- colnames(plotUMAP)[grep("res",colnames(plotUMAP))]
colData(se)$Clusters <- paste0("Cluster",as.integer(plotUMAP[,clustCol]) + 1)
colData(se)$UMAP1 <- plotUMAP$x
colData(se)$UMAP2 <- plotUMAP$y

pdf("results/LSI-Clustering-Peaks.pdf")
ggplot(plotUMAP, aes(x=x,y=y,color=res.0.8)) + geom_point(size = 0.5) + 
  theme_bw() + xlab("UMAP1") + ylab("UMAP2")
dev.off()

对于亚群分析(hematopoiesis: CD34 + bone marrow 和 DCs; tumor: T cells), 分析者重新计算了聚类的和以及log标准化的CPM值，分别从CD34+ cells 和 T cells中鉴定出前10,00和5,00个波动peak。这些波动peak用于提取稀疏二项化可及矩阵，然后做TF-IDF转换。对TF-IDF转换结果进行奇异值分解(SVD),得到更低的维度, 即用前25个维度表征原来的数据集。分贝用1-25和2-25个维度作为Seurat的输入，进行SNN graph clustering（默认resolution=0.8). 最后的结果用Seurat::RunUMAP函数进行可视化展示

细胞类群定义

以human hematopoiesis为例，介绍作者定义细胞类群的三种方法。

(1) 顺式作用元件(cis-lements, 即ATAC-seq peaks)的染色质可及性

(2) 利用和单个基因启动子相关的几个增强子的可及性计算基因活跃得分(GA)

(3) TF活跃度，通过计算每个细胞中全基因组范围内的TF结合位点可及性得出

这三种方法最大的特点就是不需要你再去测一个单细胞RNA-seq或者要一个已有的Bulk ATAC-seq谱，仅仅用当前的单细胞ATAC-seq数据即可。

这是我在哔哩哔哩上课程的文字稿版本。我有点完美主义的倾向，所以为了避免自己因为觉得自己某个部分的说话不利索，不清楚，所以总是会写稿子。

这次课程的主题是生信入门必须学习的几个Linux操作。

为什么要学习Linux

如果你要学习生物信息学，那么你一定要学会使用Linux。这是因为，绝大部分的生信软件都需要运行在Linux平台，并且有些时候，你会要用到服务器去处理大量的数据，而这些服务器几乎都是Linux环境。

当然，你可以说，我可以用云平台呀，毕竟现在很多公司都开发了云平台，我们拖拖拽拽就可以搞分析了呀。是的，现在的云平台也来越成熟，使用体验也很棒，我也很喜欢。不过我还是建议大家学Linux，因为当你学会Linux之后，你会变得更加自由，不再受限于别人的平台，而且能让日常工作变得高效起来。

如何学习Linux

那么如何学习呢？通常你问别人如何学Linux的时候，对方十有八九会给你推荐鸟哥的书，这没啥毛病，因为我也是靠鸟哥的书入门的，那本书已经被我翻成了这副模样。

学习没有捷径，但是可以少走几步路。我们不需要精通Linux，成为Linux大牛，靠Linux吃饭，我们只需要掌握一些基本操作，搞得定日常基本操作就行。就像绝大部分人用Windows的时候，也不知道如何配置Windows里的环境变量，但是一点不妨碍基本使用呀。

通过我对自己平时操作的总结，我认为学Linux起码要掌握下面这三方面内容。

• 操作迁移，就是把Windows上的常见操作都在Linux里找到对应的操作
• TAB补全，不但可以让你的操作更快，而且还能让你避免大部分的报错
• 环境变量，学习Linux一定要理解环境变量，否则你就装不好软件。

不过，最重要的是，你得先有一个Linux环境。我在哔哩哔哩上录制了两个视频，教大家如何在Windows平台上配置你的Linux环境，如何在Mac上配置Linux环境。

先看Vim视频:

• https://www.bilibili.com/video/av58673256
  再配置Linux环境：

• Windows10/Linux: https://www.bilibili.com/video/av58541835

• MacOS: https://www.bilibili.com/video/av58541835

本次视频的操作虽然是在mac上录制，但是适用于大部分的Linux环境，包括Ubuntu和CentOS，如果有什么问题，欢迎在留言区提问交流。

TAB补全

在演示Linux操作之前，我一定要先强调下tab补全的重要性，它是我认为最重要的一个操作。

初学者在学习Linux的时候，最常见的状态就是，按照教程内容，一个一个字符的敲。由于刚学习，敲代码还不熟练，因此，十有八九会出现敲错的情况，那么结果就是代码运行失败。运行失败怎么办？可能也不会去看报错，然后检查自己的输入，估计就是截图或者拍照发到群里提问了。

因此，在正式开始敲代码之前，大家先在键盘上找到tab，多敲几次，感受到它的存在。

后面，我还会不断强调它，直到它刻在你的脑子里。

操作迁移

接下来，我们来讲讲这次课程的主要内容，操作迁移。

我们在Windows上最常做的操作是什么呢？我觉得，最常见的操作应该是下面这些吧

• 浏览目录
• 切换目录
• 新建文件夹/删除文件夹
• 新建文件/删除文件
• 查看文本
• 了解资源使用情况
• 数据下载
• 安装软件

那么这些操作在Linux里的对应命令是什么呢？我们先浏览一下命令，然后通过一个模拟的项目来一个个了解他们。

• 浏览目录: ls
• 切换目录: cd, pwd
• 新建文件夹/删除文件夹: mkdir / rmdir / rm -r
• 新建文件/删除文件: vim / rm
• 查看文本: less, head, tail, cat
• 了解资源使用情况: top
• 数据下载: wget / curl
• 安装软件: apt-get/yum/conda

那么这些命令应该怎么用呢？其实很简单，就是输入命令名，后面的参数，看情况加。强调一点，如果要用到参数，那么命令名和参数之间是需要空格分隔的。

我们以一个项目为例，去学习使用最常见的Linux命令。这个任务内容如下

1, 在家目录下创建文件夹，名为 abc

2, 将abc重命名为 study

3, 下载拟南芥的注释GFF文件， 解压缩

4, 使用less查看文件内容

5, 删除gff 文件

6, 删除 study文件夹

具体操作阅读原文看我的视频。大约在5分20秒处

环境变量

环境变量也是Linux学习中非常重要的知识点，不了解他甚至都不能用好软件。

不过环境变量属于哪种你不知道，你觉得很高级，一旦知道后，却发现很简单的存在。其实我们从小学或者初中开始就通过数学了解到它。比如说x + 3 =5， x+ 4=5，其中的x就是变量，就是会变的量。

在计算机里，它就用来存储其他的值一个名字，相当于一个中介。

变量有可以分为环境变量和局部变量，环境变量就是比局部变量作用更广泛的一个变量。举个例子，当我说到爱因斯坦这个名字时，这会让你想到的是一个物理学家，但是如果你家里有一只猫，它也叫爱因斯坦，他就是一个局部变量。局部变量可以覆盖环境变量。

接下来我们将会通过一个例子，通过PATH这个环境变量来理解。

• 查看环境变量PATH

• 新建一个目录存放软件

• 将目录加入环境变量PATH

• 下载软件

• 给软件添加执行权限

• 测试软件

具体操作阅读原文看我的视频。在15分钟的时候。

视频我传到了哔哩哔哩，假如你不喜欢哔哩哔哩，想要下载我的2880 x 1800（大约3G）的视频，这是百度盘地址，链接: https://pan.baidu.com/s/17yxk78_kpeL5BA-3P0i06Q 密码:7ro0。相信我，你哔哩哔哩已经够高清了，你真的不需要20分钟3G的视频。

对于入门而言，知道的越少越好，因为细节无穷无尽。你先要上手，会敲命令了。那么后期继续深入反而就只是时间问题了。

之后，我会更新一系列入门提高的视频，介绍Linux的各种细节，比如说，

• 括号的用法
• 变量操作
• shell脚本
• …

这一部分承接从零开始学CIRCOS绘制圈图(四)，对之前的布局进行调整，使结构更加适合于发表。

染色体的位置顺序

默认情况下是会显示所有的染色体，而且会先从ath1到ath5,然后从aly1到aly8, 那么如何调整顺序呢？

需要设置的参数是chromosomes_order, 按照自己的需求调整位置。

1

chromosomes_order = aly8,aly7,aly6,aly5,aly4,aly3,aly2,aly1

虽然circos提供了一些快捷操作，比如说-,$,^, 但是远不如直接输入直观。

反转染色体的坐标顺序

默认都是从0到最大值，可以通过正则匹配的方式，将aly的染色体改成从最大值到0

1

chromosomes_reverse = /aly/

调整染色体缩放

为了让两个物种的基因组能够占Circos的两边，需要设置缩放。

1

chromosomes_scale = /aly/=0.5rn,/ath/=0.5rn

/aly/=0.5rn表示aly的所有染色体总共占据50%的空间

定义染色体的间隔

<spcaing>控制的是不同染色体之间间隔的空隙大小，比如说下面的参数表示所有染色体之间的间隔是0.005r。

1
2
3

<spacing>
default = 0.005r
</spacing>

我们设置某些染色体之间的距离大一些，比如说两个物种的起始染色体和两个物种的终止染色体

1
2
3
4
5
6
7
8
9
10
11
12

<spacing>
default = 1u

  <pairwise aly1 ath1>
   spacing = 10u
  </pairwise>
  
  <pairwise aly8 ath5>
   spacing = 10u
  </pairwise>
  
</spacing>

调整起始位置

默认情况下，ath1从-90度位置开始，由于设置了aly1和ath1之间的空隙为10u，于是看起来就不对称了。为了让结果图能够对错，需要调整角度。

1

angle_offset* = -85

最终效果

通常circos的中间部分不是空白区域，会用一条条线进行连接，表示两个染色体部分区域有关系。

数据格式

对于link，circos要求输入数据至少有6列，分别是chr1 start1 end1 chr2 start2 end2 [options]

举个例子

1

chr1	1000000	2000000	chr5	3000000	4000000

构建输入

这次会以A. lyrata 和 A.thalina的基因组为例，利用JCVI来构建Circos的输入.

新建一个项目文件夹

1

mkdir -p ath_aly && cd ath_aly

如果没有安装jcvi，用conda进行安装

1
2
3
4
5
6

conda create -n jcvi jcvi
conda activate jcvi
conda install last scipy
wget https://raw.githubusercontent.com/tanghaibao/jcvi-bin/master/bin/scip
chmod 755 scip
mv scip ~/miniconda3/bin

数据下载

从http://plants.ensembl.org/index.html分别下载这两个物种的GFF文件和CDS序列，

1
2
3
4
5
6

# Athaliana
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/fasta/arabidopsis_thaliana/cds/Arabidopsis_thaliana.TAIR10.cds.all.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.44.gff3.gz
# Alyrata
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/fasta/arabidopsis_lyrata/cds/Arabidopsis_lyrata.v.1.0.cds.all.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/gff3/arabidopsis_lyrata/Arabidopsis_lyrata.v.1.0.44.gff3.gz

数据预处理

将GFF3转成BED格式

1
2

python -m jcvi.formats.gff bed --type=mRNA --key=transcript_id Arabidopsis_thaliana.TAIR10.44.gff3.gz > ath.bed
python -m jcvi.formats.gff bed --type=mRNA --key=transcript_id Arabidopsis_lyrata.v.1.0.44.gff3.gz  > aly.bed

将BED去重复

1
2

python -m jcvi.formats.bed uniq ath.bed
python -m jcvi.formats.bed uniq aly.bed

提取CDS序列

1
2

seqkit grep -f <(cut -f 4 ath.bed ) Athaliana_167_TAIR10.cds.fa.gz | seqkit seq -i > ath.cds
seqkit grep -f <(cut -f 4 aly.bed )  Arabidopsis_lyrata.v.1.0.cds.all.fa.gz | seqkit seq -i > aly.cds

karyotype

从ENSEMBLE下载的GFF文件中，已经包含了每个基因组的大小, 当然也可以各种工具从基因组序列序列中获取大小。

karyotype.aly.txt

1
2
3
4
5
6
7
8

chr	-	aly1	aly1	0	33132539	rdylbu-11-div-1
chr	-	aly2	aly2	0	19320864	rdylbu-11-div-2
chr	-	aly3	aly3	0	24464547	rdylbu-11-div-3
chr	-	aly4	aly4	0	23328337	rdylbu-11-div-4
chr	-	aly5	aly5	0	21221946	rdylbu-11-div-5
chr	-	aly6	aly6	0	25113588	rdylbu-11-div-6
chr	-	aly7	aly7	0	24649197	rdylbu-11-div-7
chr	-	aly8	aly8	0	22951293	rdylbu-11-div-8

karyotype.tair10.txt

1
2
3
4
5

chr	-	ath1	ath1	0	30427671	brbg-10-div-1
chr	-	ath2	ath2	0	19698289	brbg-10-div-3
chr	-	ath3 	ath3	0	23459830	brbg-10-div-5
chr	-	ath4	ath4	0	18585056	brbg-10-div-7
chr	-	ath5 	ath5	0	26975502	brbg-10-div-9

因为ENSMEBLE上Athalina和Alyrata的染色体命名都是1,2,3…，就会导致CIRCOS无法正确的区分来源，因此在原本的命名前加上了物种名缩写做为标签。

同时。我们要修改之前的bed文件

1
2

sed -e 's/^/ath/' ath.uniq.bed  > ath.bed
sed -e 's/^/aly/' aly.uniq.bed  > aly.bed

links

为了构建links文件，需要利用JCVI进行共线性分析.

确保有4个文件

1
2

$ ls ???.???
aly.bed aly.cds ath.bed ath.cds

用jcvi进行分析

1
2

python -m jcvi.compara.catalog ortholog --no_strip_names ath aly
python -m jcvi.compara.synteny screen --minspan=30 --simple ath.aly.anchors ath.aly.anchors.new

其中ath.aly.anchors.simple是我们后续要用到的文件。

1
2

$ head -n 1 ath.aly.anchors.simple
AT1G24260.1	AT1G27280.1	fgenesh1_pm.C_scaffold_1002045	fgenesh2_kg.1__2877__AT1G24260.1	225	-

我们需要将ath.aly.anchors.simple里的基因名替换成实际的位置信息，将其变成符合Circos的输入信息。

我写了一个simple2links.py脚本，代码在我的GitHub上，https://github.com/xuzhougeng/myscripts。

1

python ~/myscripts/simple2links.py ath.aly.anchors.simple

最终会输出ath.aly.anchors.simple_link.txt

配置circos

接下来做如下配置。 新建一个etc文件夹，在里面建立一个links.conf用来配置links, 建立一个ticks.conf配置ticks

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

mkdir etc
# vim etc/links.conf
<links>

<link>
file          = ath.aly.anchors.simple_link.txt
radius        = 0.61r
color         = blue_a4
ribbon = yes
</link>

</links>

# vim etc/ticks.conf

show_ticks          = yes
show_tick_labels    = yes

<ticks>

radius           = 1r
color            = black
thickness        = 2p
multiplier       = 1e-6
format           = %d

<tick>
spacing        = 1u
size           = 5p
</tick>

<tick>
thickness      = 4p
spacing        = 5u
size           = 10p
show_label     = yes
label_size     = 10p
label_offset   = 10p
format         = %d
</tick>

</ticks>

编辑circos.conf

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

karyotype = karyotype.tair10.txt,karyotype.aly.txt

chromosomes_color = chr1=rdylbu-11-div-1,chr2=rdylbu-11-div-3,chr3=rdylbu-11-div-5,chr4=rdylbu-11-div-7,chr5=rdylbu-11-div-9

chromosomes_units = 1000000
<<include ./etc/ticks.conf>>

<ideogram>
<spacing>
default = 0.005r
</spacing>
radius           = 0.90r
thickness        = 20p
fill             = yes
stroke_color     = dgrey
stroke_thickness = 2p

show_label     = yes #展示label
label_font     = default # 字体
label_radius   = dims(ideogram,radius) + 0.08r #位置
label_size     = 16 # 字体大小
label_parallel = yes # 是否平行

label_format   = eval(sprintf("%s",var(chr))) # 格式
</ideogram>

<<include ./etc/links.conf>>

<image>
dir*    = .    # 输出文件夹
radius* = 500p # 图片半径
svg*    = no   # 是否输出svg
<<include etc/image.conf>>
</image>

<<include etc/colors_fonts_patterns.conf>>
<<include etc/housekeeping.conf>>

运行circos -conf circos.conf的效果如下

不难发现一个问题，里面线的颜色一模一样，不容易进行区分。虽然可以在输入ath.aly.anchors.simple_link.txt里增加颜色，但是circos提供了一个动态规则，可以更加方便的直接在配置文件里修改。

建立规则(rules)

circos的配置格式为

1
2
3
4
5
6
7
8
9
10
11
12

<rules>

<rule>
...
</rule>

<rule>
...
</rule>
...

</rules>

circos的规则可以很复杂，但是最简单的情况就是下面这种

1
2
3
4

<rule>
condition = var(chr1) eq "ath1"
color=rdylgn-5-div-1
</rule>

condition = var(chr1) eq "ath1"表示，判断link文件中左侧染色体的名字(var(chr1))是不是(eq)”ath1”，如果是的话，那么颜色就是rdylgn-5-div-1

我们可以在etc/links.conf中增加五个条件，修改后的links.conf如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

<link>
file          = ath.aly.anchors.simple_link.txt
radius        = 0.61r
color         = blue_a4
ribbon = yes

<rules>
<rule>
condition = var(chr1) eq "ath1"
color=rdylgn-5-div-1
</rule>
<rule>
condition = var(chr1) eq "ath2"
color=rdylgn-5-div-2
</rule>
<rule>
condition = var(chr1) eq "ath3"
color=rdylgn-5-div-3
</rule>
<rule>
condition = var(chr1) eq "ath4"
color=rdylgn-5-div-4
</rule>
<rule>
condition = var(chr1) eq "ath5"
color=rdylgn-5-div-5
</rule>
</rules>

</link>
</links>

运行circos -conf circos.conf的效果如下

这张图还有一个问题是，通常别人的Circos染色体都是对称排列，右边第一个是ath1，那么左边第一个也最好是aly1，那么应该如何配置呢？以及两套基因组会分别占据两侧，中间有一个明显的空隙，如何做到的呢？

这一篇会在之前的基础上开始在circos绘制基因密度信息。为了保证一致性，可以新建如下几个文件

circos.conf:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

karyotype = karyotype.tair10.txt

chromosomes_color = chr1=rdylbu-11-div-1,chr2=rdylbu-11-div-3,chr3=rdylbu-11-div-5,chr4=rdylbu-11-div-7,chr5=rdylbu-11-div-9

chromosomes_units = 1000000
<<include ticks.conf>>

<ideogram>
<spacing>
default = 0.005r
</spacing>
radius           = 0.90r
thickness        = 20p
fill             = yes
stroke_color     = dgrey
stroke_thickness = 2p

show_label     = yes #展示label
label_font     = default # 字体
label_radius   = dims(ideogram,radius) + 0.05r #位置
label_size     = 16 # 字体大小
label_parallel = yes # 是否平行

label_format   = eval(sprintf("%s",var(chr))) # 格式
</ideogram>

<image>
dir*    = .    # 输出文件夹
radius* = 500p # 图片半径
svg*    = no   # 是否输出svg
<<include etc/image.conf>>
</image>

<<include etc/colors_fonts_patterns.conf>>
<<include etc/housekeeping.conf>>

karyotype.tair10.txt

1
2
3
4
5

chr - chr1 chr1 0 30427617 chr1
chr - chr2 chr2 0 19698289 chr2
chr - chr3 chr3 0 23459830 chr3
chr - chr4 chr4 0 18585056 chr4
chr - chr5 chr5 0 26975502 chr5

ticks.conf

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

show_ticks          = yes
show_tick_labels    = yes

<ticks>

radius           = 1r
color            = black
thickness        = 2p
multiplier       = 1e-6
format           = %d

<tick>
spacing        = 1u
size           = 5p
</tick>

<tick>
thickness      = 4p
spacing        = 5u
size           = 10p
show_label     = yes
label_size     = 10p
label_offset   = 10p
format         = %d
</tick>

</ticks>

在开始之前，请确保已经安装了bedtools，如果没有的话，用conda安装

1
2

# 安装bedtools
conda install -c bioconda bedtools

数据格式

为了能够在circos绘制基因密度信息，需要先知道circos要求的输入数据格式是什么。

对于折线图(line)，散点图(scatter)，柱状图(histogram)和热图(heatmap)，Circos要求的数据输入格式相同，也就是chr start end value [options], 如果熟悉BED格式定义的话，你就会发现除了可选(options)外，Circos要求的格式就是4列的BED。

对于可选列，可以和circos.conf里的<rule>搭配使用，属于比较高级的用法。

数据预处理

为了能够获得所需的基因密度信息，我们需要下载拟南芥的GFF文件。

1
2

# download
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.44.gff3.gz

提取基因的位置信息

1

zgrep '[[:blank:]]gene[[:blank:]]' Arabidopsis_thaliana.TAIR10.44.gff3.gz | cut -f 1,4,5 | awk '{print "chr"$1"\t"$2"\t"$3}' > genes.bed

接着用bedtools以500kb为滑窗，沿染色体创建窗口

1
2

cut -d ' ' -f 3,6 karyotype.tair10.txt | tr ' ' '\t' > tair10.genome
bedtools makewindows -g tair10.genome -w 500000 > tair10.windows

最后统计信息

1

bedtools coverage -a tair10.windows -b genes.bed | cut -f 1-4 > genes_num.txt

展示信息

我们可以先用最少的参数，同时展示不同的图形。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

...
<plots>

<plot>
type      = line
thickness = 2
max_gap = 1u
file = genes_num.txt
color   = redv
r0      = 0.51r
r1      = 0.60r
</plot>

<plot>
type    = heatmap
file    = genes_num.txt
color   = spectral-5-div
r1      = 0.70r
r0      = 0.61r
</plot>

<plot>
type = scatter
fill_color       = grey
stroke_color     = black
glyph            = circle
glyph_size       = 10
file = genes_num.txt
r1   = 0.80r
r0   = 0.71r
</plot>

<plot>
type = histogram
file = genes_num.txt
r1   = 0.89r
r0   = 0.81r
</plot>

</plots>
...

运行之后，效果如下

这张图展示了巨大的进步空间，至少可以从如下几个角度进行调整，

• 柱状图有点丑，需要填充颜色
• 热图的颜色不符合直觉，基因数越少反而颜色越深
• 散点图的点太大
• 有些图形可能需要加上背景色。
• 对于一些过大的数值，最好用一种颜色表示。

当知道自己的目标后，后续的事情就是找对应参数和调整参数. 部分我用来调整参数如下

• show - 是否展示该图形
• type - 展示的图形类型
• file - 输入的数据文件所在路径
• min/max - 数据范围
• r0/r1 - 内径和外径，在圈图中的位置
• glyph - 对于散点图而言，还可以选择符号的类型，是circle, rectangle, 还是 triangle
• glyph_size - 符号的大小，单位为p
• color 散点图符号颜色，柱状图外部线的颜色
• stroke_color - 对于散点图，符号外部是否也要颜色
• stroke_thickness - 对于散点图，符号外部线的厚度
• fill_color:柱状图填充色

调整后的参数为

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

...
<plots>

<plot>
type    = heatmap
file    = genes_num.txt
color   = blues-9-seq
r1      = 0.70r
r0      = 0.61r
</plot>

<plot>
type = scatter
fill_color       = black # 填充色
stroke_color     = black
glyph            = circle
glyph_size       = 5 # 元素大小
file = genes_num.txt
r1   = 0.80r
r0   = 0.71r
</plot>

<plot>
type = histogram
file = genes_num.txt
fill_color = blue # 填充色
r1   = 0.89r
r0   = 0.81r
</plot>

</plots>
...

图形效果如下

这个效果我个人还是比较满意的。更加复杂的设置，目前还不适合我，需要一步一步来。

下一部分将会介绍如何在circos展示两个物种的共线性，也就是circos不同染色体之间的连线。

在从零开始学CIRCOS绘制圈图(一)， 我们已经绘制出一个比较丑的circos图，这一部分是讲解一些细节。

这一部分会从上一步的两个文件开始，分别是

karyotype.tair10.txt

1
2
3
4
5

chr - chr1 chr1 0 30427617 chr1
chr - chr2 chr2 0 19698289 chr2
chr - chr3 chr3 0 23459830 chr3
chr - chr4 chr4 0 18585056 chr4
chr - chr5 chr5 0 26975502 chr5

circos.conf

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

karyotype = karyotype.tair10.txt

<ideogram>
<spacing>
default = 0.005r
</spacing>
radius           = 0.90r
thickness        = 20p
fill             = yes
stroke_color     = dgrey
stroke_thickness = 2p
</ideogram>

<image>
<<include etc/image.conf>>
</image>

<<include etc/colors_fonts_patterns.conf>>
<<include etc/housekeeping.conf>>

配色

第一个要解决的问题是配色问题，其实上一部分的配色定义依旧有点丑。如果不想重新打开文件修改的话，其实除了在karyotype文件里定义颜色外，我们还可以直接在circos.conf文件定义颜色.

举个例子，在karyotype = karyotype.tair10.txt后加一行

1

chromosomes_color = chr1=rdylbu-11-div-1,chr2=rdylbu-11-div-3,chr3=rdylbu-11-div-5,chr4=rdylbu-11-div-7,chr5=rdylbu-11-div-9

问题是，我们怎们知道这些名字后所代表的颜色呢？

Circos中颜色的命名格式为PALETTE-NUMCOLORS-TYPE-IDX:

• PALETTE:调色版名，如rdylbu
• NUMCOLORS: 颜色数目, 11
• 调色版类型: div(diverging), seq(sequential), qual(qualitative)
• IDX: 调色版中的颜色索引

而Circos颜色来自于http://colorbrewer2.org

因此，gnbu-9-seq对应的是就是下图的9-class GnBu

配色不仅仅用在karyotype中，在后续的热图和柱状图等还会涉及到它，毕竟一张好看的图，配色占了很大的比例。

显示标签

默认输出图片是没有染色体名字的标签，需要在<ideogram>里添加和label有关的参数

1
2
3
4
5
6
7
8
9
10
11

...
</spacing>
...
show_label     = yes #展示label
label_font     = default # 字体
label_radius   = dims(ideogram,radius) + 0.05r #位置
label_size     = 16 # 字体大小
label_parallel = yes # 是否平行

label_format   = eval(sprintf("%s",var(chr))) # 格式
</ideogram>

关于标签(label)， 更详细的介绍在http://circos.ca/documentation/tutorials/ideograms/labels/

重新运行之后，发现字相对而言太小了。有两种结局方案，一种是调整label_size，比如说48，另一种是调整图片整体大小。

输出设置

默认情况下，输出图片的 半径是1500p, 所以设置的label_size=16就会显得特别小。我们可以在配置文件中调整图片的班级，以及其他设置。

1
2
3
4
5
6

<image>
dir*    = .    # 输出文件夹
radius* = 500p # 图片半径
svg*    = no   # 是否输出svg
<<include etc/image.conf>>
</image>

注: 在参数后加一个*表示覆盖已有的设置，比如说svg*=no就是覆盖已有的svg=yes。

运行结果如下

刻度(ticks)

大部分的circos还会有刻度来展示染色体的大小。由于ticks的定制比较复杂，所以一般会单独搞一个配置文件，ticks.conf存放ticks相关的参数设置.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

chromosomes_units = 1000000
show_ticks          = yes
show_tick_labels    = yes

<ticks>

radius           = 1r
color            = black
thickness        = 2p

multiplier       = 1e-6 #输出的标签为实际长度与其相乘

format           = %d # %d表示显示整数

<tick>
spacing        = 1u
size           = 5p
</tick>

<tick>
spacing        = 5u
size           = 10p
show_label     = yes
label_size     = 10p
label_offset   = 10p
format         = %d
</tick>

</ticks>

用 show_ticks和 show_tick_labels控制是否展示刻度，以及刻度对应的标签。

<ticks与</ticks>里控制刻度的全局参数，例如位置为1r(radius=1r), 颜色为黑色(color=black), 厚度为2p(thickness=2p), 由于默认直接展示染色体的实际位置，因此会显示1000000这种结果，所以定义multiplier=1e-6, 实际显示结果为 1000000 * 1e-6 = 1。

后面就可以通过<tick>和 </tick>来分别绘制不同类型的tick，重要的参数如下：

• spacing表示刻度之间的距离，1u表示一个长度单位，需要在circos.conf文件里通过chromosome_unit来定义，通常都是chromosome_unit=1000000。
• size表示tick的长度
• show_label: 控制是否展示标签，默认不展示。
• label_offset: 则是让label往外在偏移一点距离

circos -conf circos.conf运行结果如下

我们发现有些刻度的标签和染色体标签发生了重叠，这个可以通过label_radius进行调整。

单位

上面出现了控制图形不同元素大小的三个单位，p,r,u。p(pixels), 表示绝对大小， r(relative), 相对大小， u(chromosome unit)。 如果使用p作为单位，需要考虑最终输出图形<image>定义的radius。 而r是相对大小，会随着最终图形大小而发生变换。u一般在显示刻度时使用。

这一部分是在原来简陋的输出上进行了美化，没有用到除了染色体长度以外的信息。下一部分介绍如何展示基因密度等信息。

在如何从BAM文件中提取fastq文章里其实也发现了从BAM里面提取Fastq是有些麻烦，只不过最后通过samtools的子命令实现了数据提取，实现功能之后也没有再去思考如何提高效率。

最近读到每周文献-190419-植物单细胞BAM重比对以及假基因研究时，发现里面提到了一个工具叫做 bazam, 功能就是提取Fastq文件，文章发表在 Genome Biology 上。

软件地址为https://github.com/ssadedin/bazam，因为他是个Groovy工具（据说Groovy比Java好写）所以安装很方便，

1
2
3
4

git clone https://github.com/ssadedin/bazam.git
cd bazam
git submodule update --init --recursive
./gradlew clean jar

可以通过如下命令检查是否安装成功

1

java -jar build/libs/bazam.jar -bam  test.bam  > tmp.fastq

假如你原来的BAM文件是双端测序，想提取成两个文件，那么可以用如下命令。注：这里的your.bam 指的是你实际BAM文件

1
2

java -Xmx16G -jar build/libs/bazam.jar \
    -bam your.bam -r1 your_R1.fq -r2 your_R2.fq &

此外，bazam还支持使用-L指定区间来提取某个区间的数据。

尽管是可以从BAM文件中提取Fastq文件，但如果原始BAM文件过滤了未必对的read，那么你还是会损失掉这部分信息，建议保存好原来的Fastq文件。

这里补充下为啥从BAM文件中提取双端序列那么麻烦，这是因为一般而言BAM文件都是按照位置信息排序，想要找到配对的reads，要么是根据read的编号进行排序（这个方法要求额外的内存和存储空间），或者就是在提取的时候记录当前的read的ID，再找到另一端ID后释放内存空间。

参考资料

• Bazam: a rapid method for read extraction and realignment of high-throughput sequencing data
• https://github.com/ssadedin/bazam

虽然高通量测序分析最常用的操作是将fastq比对到参考基因组得到BAM文件，但偶尔我们也需要提取BAM文件中特定区域中fastq。最开始我认为这是一个非常简单的操作，因为samtools其实已经提供了相应的工具samtools fastq.

以biostar handbook的Ebola病毒数据为例，首先获取比对得到的BAM文件。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

# 建立文件夹
mkdir -p refs
# 根据Accession下载
ACC=AF086833
REF=refs/$ACC.fa
efetch -db=nuccore -format=fasta -id=$ACC | seqret -filter -sid $ACC > $REF
bwa index $REF
SRR=SRR1553500
fastq-dump -X 100000 --split-files $SRR
BAM=$SRR.bam
R1=${SRR}_1.fastq
R2=${SRR}_2.fastq
TAG="@RG\tID:$SRR\tSM:$SRR\tLB:$SRR"
bwa mem -R $TAG $REF $R1 $R2 | samtools sort > $BAM
samtools index $BAM

让我们尝试提取下前1kb比对的fastq文件

1
2
3

samtools view -h SRR1553500.bam KJ660346.1:1-1000 | samtools fastq -1 read_1.fq -2 read_2.fq -
#[M::bam2fq_mainloop] discarded 0 singletons
#[M::bam2fq_mainloop] processed 8702 reads

看起来似乎我们成功了，那让我们把这些序列比对回去吧

1
2

bwa mem refs/AF086833.fa read_1.fq read_2.fq
[mem_sam_pe] paired reads have different names: "SRR1553500.11772", "SRR1553500.43775"

你会发现这里产生了错误。这个错误说的是配对的read名字不同。导致问题的原因是：samtools fastq会按照fastq在bam里原本的顺序进行提取，而按照位置排序的bam的两个配对read并不是前后紧挨着的。所以我们需要进行一波排序。

1
2
3
4
5
6
7
8
9
10

cat read_1.fq \
| paste - - - - \
| sort -k1,1 -S 3G \
| tr '\t' '\n' \
| gzip > read_1.fq.gz
cat read_2.fq \
| paste - - - - \
| sort -k1,1 -S 3G \
| tr '\t' '\n' \
| gzip > read_2.fq.gz

这里将5个shell命令用管道连起来完成这个任务，咋看起来比较复杂，但是理解每一步的目的，每个参数的作用，以后就能活用起来了。

首先用cat逐行读取read，传入到paste中。paste的作用是将多个文件合并成一个文件，做法比较简单粗暴，举个例子

1
2
3
4
5
6
7
8
9
10
11
12

cat A.txt
A
B
C
cat B.txt
D
E
F
paste A.txt B.txt
A	D
B	E
C	F

cat read_1.fq | paste - - - - 表示读取4行合并成一行。相当于fastq里的read都会用一行而不是四行表示。接着用sort按照第一行(-k1,1)排序，内存的缓冲数据为3G(-S 3G). 然后将制表符替换成换行符，让一行又变回到4行，后续是压缩成gz，降低磁盘占用。

这样做其实还没又解决问题，因为这两个文件里的reads数其实不相同，我们需要进一步从里面剔除掉只有一方才有的read。这里使用的seqkit，一个非常强大的序列处理轮子。

1
2
3
4
5
6

seqkit grep -f <(seqkit seq -n -i read_1.fq.gz ) read_2.fq.gz  > read_flt_2.fq
seqkit grep -f <(seqkit seq -n -i read_2.fq.gz ) read_1.fq.gz  > read_flt_1.fq
wc -l read_flt_1.fq
   15028 read_flt_1.fq
wc -l read_flt_2.fq
   15028 read_flt_2.fq

终于，你得到了配对存放的两个fastq文件。

通过这个案例我想表达的观点是：虽然有些时候没有现成的工具可以解决你的问题，但是你可以通过组合几个现有的工具来完成。也就是将一个大问题拆解成几个独立的小问题，然后逐个击破。

当然，你首先得熟悉一些已有的工具，比如说常用unix命令:sort, paste, tr等，和一些非常强大生信命令行工具：samtools, bedtools, seqkit等。

一个更佳方便的方法

前面只是问题的一种方法，可能还不是最好的，但是我当初最先想到的一个比较粗暴的方法。当然更佳巧妙的方法是，提取序列之后可以按照read name排序，然后提取

1

samtools view -h SRR1553500.bam KJ660346.1:1-1000 | samtools sort -n | samtools fastq -1 read_1.fq -2 read_2.fq -s singleton.fq -

是不是更好理解了。