Allele.ctg.table是ALLHiC用于过滤多倍体基因组中因等位序列相似引起的HiC噪音的必要输入。

构建的方法有很多种，这里列举两个方法。

方法一: 基于BLAST

基于BLAST的方法需要你提供一个染色体级别组装的近缘物种。需要提供4个文件，近缘物种的CDS和GFF文件，多倍体物种的CDS和GFF文件。最重要的是，你得保证CDS序列中的编号和GFF文件的基因编号一致，这也是后续处理至关重要的一步。

问题一: 我们如何快速获取待组装物种的CDS和GFF文件？最快速的方法就是利用转录组序列基于STAR + StringTie的组合进行拼接，并利用gffread从中提取出cdna序列。

假设你的多倍体文件名为contig.fa, 转录组数据为read_R1.fq.gz, read_R2.fq.gz

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

mkdir -p STAR
ref=contig.fa
THREADS=40
fq1=read_R1.fq.gz
fq2=read_R2.fq.gz
prefix=sample
# 建立索引
STAR \
  --runThreadN $THREADS \
  --runMode genomeGenerate \
  --genomeDir STAR \
  --genomeFastaFiles $ref
#比对
STAR \
	--genomeDir STAR \
	--runThreadN $THREADS \
	--readFilesIn $fq1 $fq2 \
	--readFilesCommand zcat \
	--outFileNamePrefix $prefix\
	--outSAMtype BAM SortedByCoordinate \
	--outBAMsortingThreadN $THREADS \
	--outSAMstrandField intronMotif \
	--outFilterIntronMotifs RemoveNoncanonical
# 拼接转录本
stringtie ${prefix}Aligned.sortedByCoord.out.bam -p $THREADS -o qry.gtf
# GTF 转成 GFF3
gffread qry.gtf -o qry.gff
# 提取cdna
gffread qry.gff -g $ref -w qry.fa

问题二: 后续的处理脚本要求输入序列的编号和GFF的基因编号相同，并且编号在GFF文件中的格式为NAME=xxx。

对于我们上一步处理的qry.gff和qry.fa，我们只需要简单粗暴地替换qry.gff里的信息即可。

1

sed -e 's/transcript/gene/' -e 's/ID/Name/' qry.gff > qry_gene.gff

对于近缘物种的cds序列和gff文件，考虑到后续的脚本只认gene所在的行，而cds里面存放的是转录本序列，我们可以提取mRNA所在行，然后将其改名成gene。假设参考序列是ref.fa, 参考的GFF文件名是ref.gff

1
2

ref=ref.gff
grep '[[:blank:]]mRNA[[:blank:]]' $ref | sed -e 's/mRNA/gene/' -e 's/ID/Name/' > ref_gene.gff

综上你得到了如下四个文件, ref.fa, ref_gene.gff, qry.fa, qry_gene.gff

下面的操作，参考自官方教程完成，

第一步： 建立BLAST索引

1

makeblastdb -in ref.fa -dbtype nucl

第二步: 将我们物种序列比对到近缘物种

1

blastn -query qry.fa -db ref.fa -out qry_vs_ref.blast.out -evalue 0.001 -outfmt 6 -num_threads 4 -num_alignments 1

第三步: 过滤低质量的HIT，此处标准是相似度低于60%, 覆盖度低于80%

1

blastn_parse.pl -i qry_vs_ref.blast.out -o Eblast.out -q qry.fa -b 1 -c 0.6 -d 0.8 

第四步: 基于BLAST结果鉴定等位序列

1
2

classify.pl -i Eblast.out -p 4 -r ref_gene.gff -g qry_gene.gff
# -p 是等位基因的数目

最终输出的Allele.ctg.table就是后续ALLHiC的prune输入。

方法二: 基于GMAP

第二个方法比较简单，只需要提供多倍体基因组的基因组序列和近缘物种的cds序列即可

第一步: 运行GMAP获取GFF3文件

1
2
3
4
5

qry=target.genome  # 你的多倍体序列
ref_cds=reference.cds.fasta # 同源物种的cds
N=4 #你的基因组倍性, 4表示的是4倍体
gmap_build -D . -d DB $qry
gmap -D . -d DB -t 12 -f 2 -n $N  $ref_cds > gmap.gff3

第二步: 获取 allelic.ctg.table

1

perl gmap2AlleleTable.pl referenece.gff3

gmap2AlleleTable.pl脚本是ALLHiC的一部分。

参考资料

• https://github.com/tangerzhang/ALLHiC/issues/16
• https://github.com/tangerzhang/ALLHiC/wiki/ALLHiC:-identify-allelic-contigs

基因组组装大致可以分为三步（1）根据序列之间的重叠情况构建出contig，（2）基于二代的mate pair文库或光学图谱将contig搭建成scaffold，（3）对scaffold进行排序和调整方向得到最终的准染色体级别的基因组。

目前的三代测序组装能够搞定第一步和第二步。而在将contig/scaffold提升至准染色体水平上，有4种方案可选。一种是基于遗传图谱，一种是利用BioNano DLS光学图谱，一种是利用近缘物种的染色体同源性，还有一种就是HiC。其中HiC技术是三者中较为简单的一个，不需要高质量的DNA文库，也不需要一个很大的群体，结果也比较准确可信。

HiC的文库构建示意图如下，我们所需要的就是最终双端测序的两端序列之间的距离关系。

目前利用HiC数据进行组装软件有 LACHESIS, HiRise, SALSA1, 3D-DNA等，这些软件在动物基因组上和简单植物基因组上表现都不错，但是不太适合直接用于多倍体物种和高杂合物种的组装上。主要原因就是等位基因序列的相似性，使得不同套染色体之间的contig出现了假信号，最终错误地将不同套染色体的contig连在了一起。最近在Nature Plants发表的ALLHiC流程就是用来解决多倍体物种和高杂合度基因组的HiC组装难题。

ALLHiC流程一览

ALLHiC一共分为五步(见下图，Zhang et al., 2019)，pruning, partition, rescue, optimization, building，要求的输入文件为HiC数据比对后的BAM和一个Allele.ctg.table。

其中pruning步骤是ALLHiC区别于其他软件的关键一步。因此我专门将其挑选出来进行介绍，红色实线是潜在的坍缩区域（组装时因为序列高度相似而没有拆分），而其他颜色实线则是不同的单倍型（我用浅灰色椭圆进行区分）。粉红色虚线指的是等位基因间的HiC信号，而黑色虚线则是坍缩区域和未坍缩区域的HiC信号。

ALLHiC在这一步会根据提供的Allele.ctg.table过滤BAM文件中等位基因间的HiC信号，同时筛选出坍缩区域和未坍缩区域的HiC信号。这些信号会用于Rescue步骤，将未锚定contig分配到已分组的contigs群。

软件安装

ALLHiC的安装非常简单，按照习惯，我将软件安装在~/opt/biosoft下

1
2
3
4
5
6
7
8
9

mkdir -p ~/opt/biosoft && cd ~/opt/biosoft 
git clone https://github.com/tangerzhang/ALLHiC
cd ALLHiC
mv allhic.v0.9.8 bin/allhic
chmod +x bin/*
chmod +x scripts/*  
# 添加到环境变量
PATH=$HOME/opt/biosoft/ALLHiC/scripts/:$HOME/opt/biosoft/ALLHiC/bin/:$PATH
export PATH

此外ALLHiC还依赖于samtools(v1.9), bedtools 和 Python 3环境的matplotlib(v2.0+)，这些可以通过conda一步搞定。

1

conda create -y -n allhic python=3.7 samtools bedtools matplotlib

之后检查下是否成功安装

1
2
3

$ conda activate allhic
$ allhic -v
$ ALLHiC_prune

你可能会遇到如下的报错

1

ALLHiC_prune: /lib64/libstdc++.so.6: version `GLIBCXX_3.4.21' not found (required by ALLHiC_prune)

这是GLIBC过低导致，但是不要尝试动手去升级GLIBC（你承担不起后果的），conda提供了一个比较新的动态库，因此可以通过如下方法来解决问题

1

export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:~/opt/miniconda3/lib

ALLHiC分析实战

感谢张兴坦老师提供的测试数据，我将其中的BAM转成了原始FastQ格式，以便从头讲解。数据在百度网盘,提取码：17yk

输入文件需要有4个，contig序列，等位contig信息和两个双端测序数据

1
2
3
4
5

$ ls -1
Allele.ctg.table
draft.asm.fasta
reads_R1.fastq.gz
reads_R2.fastq.gz

第一步: 建立索引

1
2

samtools faidx draft.asm.fasta 
bwa index -a bwtsw draft.asm.fasta  

第二步: 序列回贴。这一步的限速步骤是bwa sampe，因为它没有多线程参数。如果数据量很大，可以先将原始的fastq数据进行拆分，分别比对后分别执行bwa sampe，最后合并成单个文件。

1
2
3

bwa aln -t 24 draft.asm.fasta reads_R1.fastq.gz > reads_R1.sai  
bwa aln -t 24 draft.asm.fasta reads_R2.fastq.gz > reads_R2.sai  
bwa sampe draft.asm.fasta reads_R1.sai reads_R2.sai reads_R1.fastq.gz reads_R2.fastq.gz > sample.bwa_aln.sam  

第三步: SAM预处理，移除冗余和低质量信号，提高处理效率

1
2
3
4
5

PreprocessSAMs.pl sample.bwa_aln.sam draft.asm.fasta MBOI
# 如果已有BAM文件
# PreprocessSAMs.pl sample.bwa_aln.bam draft.asm.fasta MBOI
filterBAM_forHiC.pl sample.bwa_aln.REduced.paired_only.bam sample.clean.sam
samtools view -bt draft.asm.fasta.fai sample.clean.sam > sample.clean.bam

其中filterBAM_forHiC.pl的过滤标准是比对质量高于30(MQ), 只保留唯一比对(XT:A:U), 编辑距离(NM)低于5, 错误匹配低于(XM)4, 不能有超过2个的gap(XO,XG)

第四步(可选): 对于多倍体或者是高杂合的基因组，因为等位基因的序列相似性高，那么很有可能会在不同套基因组间出现假信号，因此需要构建Allele.ctg.table, 用于过滤这种假信号。

1

ALLHiC_prune -i Allele.ctg.table -b sample.clean.bam -r draft.asm.fasta  

这一步生成prunning.bam用于后续分析

第五步: 这一步是根据HiC信号将不同的contig进行分组，分组数目由-k控制。如果跳过了第四步，那么可以直接用第三步的结果sample.clean.bam

1

ALLHiC_partition -b prunning.bam -r draft.asm.fasta -e AAGCTT -k 16  

这一步会生成一系列以prunning开头的文件

• 分组信息: prunning.clusters.txt
• 每个分组对应的contig: prunning.counts_AAGCTT.XXgYY.txt:
• 每个contig长度和count数: prunning.counts_AAGCTT.txt

第六步: 将未锚定的contig分配已有的分组中。

1
2
3

ALLHiC_rescue -b sample.clean.bam -r draft.asm.fasta \
    -c prunning.clusters.txt \
    -i prunning.counts_AAGCTT.txt

这一步根据之前prunning.counts_AAGCTT.XXgYY.txt对应的groupYY.txt

第七步: 优化每一组中的contig的排序和方向

1
2
3
4
5
6

# 生成.clm文件
allhic extract sample.clean.bam draft.asm.fasta --RE AAGCTT  
# 优化
for i in group*.txt; do
    allhic optimize $i sample.clean.clm
done

这一步会基于groupYY.txt生成对应的groupYY.tour

第八步: 将tour格式转成fasta格式，并生成对应的agp。

1

ALLHiC_build draft.asm.fasta  

这一步生成两个文件，groups.asm.fasta和groups.agp。其中groups.asm.fasta就是我们需要的结果。

第九步: 构建染色质交互矩阵，根据热图评估结果

1
2
3

samtools faidx groups.asm.fasta
cut -f 1,2 groups.asm.fasta.fai  > chrn.list
ALLHiC_plot sample.clean.bam groups.agp chrn.list 500k pdf

使用ALLHiC的几个注意事项:

1. ALLHiC依赖于初始的contig，如果嵌合序列和坍缩序列比例过高，那么ALLHiC结果也会不准确。根据文章，ALLHiC能够处理~10%的嵌合比例，~20%的坍缩比例。因此最好是用类似于Canu这种能够区分单倍型的组装软件。
2. 单倍型之间序列相似度不能太高，否则会出现大量非唯一比对，降低可用的HiC信号
3. 构建Allele.ctg.table需要一个比较近缘的高质量基因组
4. 不要用过短的contig，因为短的contig信号少，很容易放到错误的区域
5. K值的设置要根据实际的基因组数目设置，如果你发现输出结果中某些group过大，可以适当增大k值

参考资料

• https://github.com/tanghaibao/allhic
• https://github.com/tangerzhang/ALLHiC/wiki
• Zhang, X., Zhang, S., Zhao, Q., Ming, R., and Tang, H. (2019). Assembly of allele-aware, chromosomal-scale autopolyploid genomes based on Hi-C data. Nat. Plants 5, 833–845.
• Zhang, J., Zhang, X., Tang, H., Zhang, Q., Hua, X., Ma, X., Zhu, F., Jones, T., Zhu, X., Bowers, J., et al. (2018). Allele-defined genome of the autopolyploid sugarcane Saccharum spontaneum L. Nature Genetics 50, 1565.

NextPolish是武汉未来组开发的一个三代基因组polish工具（另外一个常用软件是Pilon）。NextPolish可以使用二代短读序列或者三代序列或者两者结合去纠正三代长读长序列在组装时导致的碱基错误(SNV/Indel)。由于它是专为polish设计，因此在运行速度和内存使用上都优与Pilon。

软件安装

先确保自己的服务器上安装了Python2.7， 且有Shutil和Signal，或者你可以利用conda新建一个python2.7的环境。

1
2
3
4
5
6

# shell
python -V
Python 2.7.15
# Python 交互命令行
import shutil
import signal

之后到https://github.com/Nextomics/NextPolish/releases找最新的版本，写这篇时的最新版本就是1.05

1
2
3
4
5
6
7
8

mkdir -p ~/opt/biosoft
cd ~/opt/biosoft
wget https://github.com/Nextomics/NextPolish/releases/download/v1.0.5/NextPolish.tgz
tar -zxvf NextPolish.tgz
# 编译软件
cd NextPolish && make -j 10
# 加入到.bashrc或.zshrc
export PATH=~/opt/biosoft/NextPolish:$PATH

软件使用

注意：如果你的基因组用的是miniasm这类缺少consensus步骤的组装软件，那么你需要先用运行如下命令，或者是运行racon利用三代序列进行polish。否则，由于基因组上存在过高的错误率，导致二代序列错误比对，影响polish效果。

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

threads=20  
genome=input.genome.fa # 组装的基因组
lgsreads=input.lgs.reads.fq.gz # 三代长度序列
# 将三代回帖到参考基因组
minimap2 -a -t ${threads} -x map-ont/map-pb ${genome} ${lgsreads}| \
    samtools view -F 0x4 -b - | \
    samtools sort - -m 2g -@ ${threads} -o genome.lgs.bam  
#建立索引
samtools index -@ ${threads} genome.lgs.bam
samtools faidx ${genome}
# 使用nextPolish.py 进行polish
python ~/opt/biosoft/NextPolish/lib/nextPolish.py \
    -g ${genome} -t 5 --bam_lgs genome.lgs.bam -p ${threads} > genome.lgspolish.fa

生成的genome.lgspolish.fa才能用于后续的二代polish步骤。

NextPolish要求我们需要准备两个文件：

• run.cfg: 配置文件，
• sgs.fofn: 二代测序文件的位置信息

以使用NextDenovo组装Nanopore数据文章组装的结果为例进行介绍。在分析目录下有三个文件。

• 三代组装结果: nextgraph.assembly.contig.fasta
• 二代序列: ERR2173372_1.fastq,ERR2173372_2.fastq

第一步：创建一个文件，用于记录二代序列的位置信息

1

realpath ERR2173372_1.fastq ERR2173372_2.fastq  > sgs.fofn

第二步：配置run.cfg文件

1
2

# 从NextPolish目录下复制配置文件
cp ~/opt/biosoft/NextPolish/doc/run.cfg run2.cfg

修改配置文件

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

[General]
job_type = local
job_prefix = nextPolish
task = default
rewrite = 1212
rerun = 3
parallel_jobs = 2
multithread_jobs = 10
genome = ./nextgraph.assembly.contig.fasta
genome_size = auto
workdir = ./01_rundir
polish_options = -p {multithread_jobs}

[sgs_option]
sgs_fofn = ./sgs.fofn
sgs_options = -max_depth 100

其中需要修改的参数为，其余参数查看官方的参数配置说明:

• job_type: 任务类型，local表示单个节点运行。由于NextPolish使用DRMAA进行任务投递，因此还支持，SGE, PBS和SLURM
• task: 任务类型， 用12,1212,121212,12121212来设置polish的轮数，建议迭代2轮就可以了。
• parallel_jobs和multithread_jobs表示同时投递的任务数和每个任务的线程数，此处2 X 10=20
• genome: 表示组装基因组的位置
• workdir: 输出文件所在目录
• sgs_options: 该选项设置二代测序polish的参数，包括-use_duplicate_reads, -unpaired, -max_depth, -bwa, -minimap2(默认使用)

运行方法

1

nextPolish run2.cfg &

在最后输出日志中，会提示最终存放的文件在什么位置，然后将这些文件合并到单个文件即可。

实用Seurat自带的热图函数DoHeatmap绘制的热图，感觉有点不上档次，于是我尝试使用ComplexHeatmap这个R包来对结果进行展示。

个人觉得好的热图有三个要素

• 聚类: 能够让别人一眼就看到模式
• 注释: 附加注释能提供更多信息
• 配色: 要符合直觉，比如说大部分都会认为红色是高表达，蓝色是低表达

在正式开始之前，我们需要先获取一下pbmc的数据，Seurat提供了R包SeuratData专门用于获取数据

1
2
3

devtools::install_github('satijalab/seurat-data')
library(SeuratData)
InstallData("pbmc3k")

加载数据并进行数据预处理，获取绘制热图所需的数据

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

library(SeuratData)
library(Seurat)
data("pbmc3k")
pbmc <- pbmc3k
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)

pbmc.markers <- FindAllMarkers(pbmc, 
                               only.pos = TRUE, 
                               min.pct = 0.25, 
                               logfc.threshold = 0.25)

先感受下Seurat自带热图

1
2

top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()

下面则是介绍如何用R包ComplexHeatmap进行组图，虽然这个R包名带着Complex，但是并不是说这个R包很复杂，这个Complex应该翻译成复合，也就是说这个R包能在热图的基础上整合很多信息。

先安装并加载R包。

1
2

BiocManager::install("ComplexHeatmap")
library(ComplexHeatmap)

为了手动绘制一个热图，要从Seurat对象中提取所需要的表达量矩阵。我提取的是原始的count值，然后用log2(count + 1)的方式进行标准化

1
2

mat <- GetAssayData(pbmc, slot = "counts")
mat <- log2(mat + 1)

获取基因和细胞聚类信息

1
2

gene_features <- top10
cluster_info <- sort(pbmc$seurat_annotations)

对表达量矩阵进行排序和筛选

1

mat <- as.matrix(mat[top10$gene, names(cluster_info)])

用Heatmap绘制热图。对于单细胞这种数据，一定要设置如下4个参数

• cluster_rows= FALSE: 不作行聚类
• cluster_columns= FALSE: 不作列聚类
• show_column_names=FALSE: 不展示列名
• show_row_names=FALSE: 不展示行名，基因数目不多时候可以考虑设置为TRUE

1
2
3
4
5

Heatmap(mat,
        cluster_rows = FALSE,
        cluster_columns = FALSE,
        show_column_names = FALSE,
        show_row_names = TRUE)

从图中，我们可以发现以下几个问题：

• 长宽比不合理，当然这和绘图函数无关，可以在保存时修改长宽比
• 基因名重叠，考虑调整大小，或者不展示，或者只展示重要的基因
• 颜色可以调整
• 缺少聚类信息

这些问题，我们可以通过在ComplexHeatmap Complete Reference查找对应信息来解决。

配色方案

在热图中会涉及到两类配色，一种用来表示表达量的连续性变化，一种则是展示聚类。有一个神奇的R包就是用于处理配色，他的Github地址为https://github.com/caleblareau/BuenColors。

1
2

devtools::install_github("caleblareau/BuenColors")
library("BuenColors")

它提供了一些列预设的颜色，比方说jdb_color_maps

1
2
3
4
5
6

      HSC       MPP      LMPP       CMP       CLP       MEP       GMP 
"#00441B" "#46A040" "#00AF99" "#FFC179" "#98D9E9" "#F6313E" "#FFA300" 
      pDC      mono     GMP-A     GMP-B     GMP-C       Ery       CD4 
"#C390D4" "#FF5A00" "#AFAFAF" "#7D7D7D" "#4B4B4B" "#8F1336" "#0081C9" 
      CD8        NK         B 
"#001588" "#490C65" "#BA7FD0"

这些颜色就能用于命名单细胞的类群，比如说我选择了前9个

1
2

col <- jdb_color_maps[1:9]
names(col) <- levels(cluster_info)

增加列聚类信息

Heatmap的row_split和column_split参数可以通过设置分类变量对热图进行分隔。更多对热图进行拆分，可以参考Heatmap split

1
2
3
4
5
6

Heatmap(mat,
        cluster_rows = FALSE,
        cluster_columns = FALSE,
        show_column_names = FALSE,
        show_row_names = FALSE,
        column_split = cluster_info)

只用文字描述可能不够好看，最好是带有颜色的分块图，其中里面的颜色和t-SNE或UMAP聚类颜色一致，才能更好的展示信息。

为了增加聚类注释，我们需要用到HeatmapAnnotation函数，它对细胞的列进行注释，而rowAnnotation函数可以对行进行注释。这两个函数能够增加各种类型的注释，包括条形图，点图，折线图，箱线图，密度图等等，这些函数的特征是anno_xxx，例如anno_block就用来绘制区块图。

1
2
3
4

top_anno <- HeatmapAnnotation(
  cluster = anno_block(gp = gpar(fill = col), # 设置填充色
                       labels = levels(cluster_info), 
                       labels_gp = gpar(cex = 0.5, col = "white"))) # 设置字体

其中anno_block中的gp参数用于设置各类图形参数，labels设置标签，labels_gp设置和标签相关的图形参数。可以用?gp来了解有哪些图形参数。

1
2
3
4
5
6
7
8

Heatmap(mat,
        cluster_rows = FALSE,
        cluster_columns = FALSE,
        show_column_names = FALSE,
        show_row_names = FALSE,
        column_split = cluster_info,
        top_annotation = top_anno, # 在热图上边增加注释
        column_title = NULL ) # 不需要列标题

突出重要基因

由于基因很多直接展示出来，根本看不清，我们可以强调几个标记基因。用到两个函数是rowAnnotation和anno_mark

已知不同类群的标记基因如下

我们需要给anno_mark提供基因所在行即可。

1
2
3
4
5
6

mark_gene <- c("IL7R","CCR7","IL7R","S100A4","CD14","LYZ","MS4A1","CD8A","FCGR3A","MS4A7","GNLY","NKG7","FCER1A", "CST3","PPBP")
gene_pos <- which(rownames(mat) %in% mark_gene)

row_anno <-  rowAnnotation(mark_gene = anno_mark(at = gene_pos, 
                                     labels = mark_gene))

接着绘制热图

1
2
3
4
5
6
7
8
9

Heatmap(mat,
        cluster_rows = FALSE,
        cluster_columns = FALSE,
        show_column_names = FALSE,
        show_row_names = FALSE,
        column_split = cluster_info,
        top_annotation = top_anno,
        right_annotation = row_anno,
        column_title = NULL)

关于如何增加标记注释，参考mark-annotation

调增图例位置

目前的热图还有一个问题，也就是表示表达量范围的图例太占位置了，有两种解决方法

• 参数设置show_heatmap_legend=FALSE直接删掉。
• 利用heatmap_legend_param参数更改样式

我们根据legends这一节的内容进行一些调整

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

Heatmap(mat,
        cluster_rows = FALSE,
        cluster_columns = FALSE,
        show_column_names = FALSE,
        show_row_names = FALSE,
        column_split = cluster_info,
        top_annotation = top_anno,
        right_annotation = row_anno,
        column_title = NULL,
        heatmap_legend_param = list(
          title = "log2(count+1)",
          title_position = "leftcenter-rot"
        ))

因为ComplextHeatmap是基于Grid图形系统，因此可以先绘制热图，然后再用grid::draw绘制图例，从而实现将条形图的位置移动到图中的任意位置。

先获取绘制热图的对象

1
2
3
4
5
6
7
8
9
10
11

p <- Heatmap(mat,
        cluster_rows = FALSE,
        cluster_columns = FALSE,
        show_column_names = FALSE,
        show_row_names = FALSE,
        column_split = cluster_info,
        top_annotation = top_anno,
        right_annotation = row_anno,
        column_title = NULL,
        show_heatmap_legend = FALSE
        )

根据p@matrix_color_mapping获取图例的颜色的设置，然后用Legend构建图例

1
2
3
4
5
6
7
8
9
10
11

col_fun  <- circlize::colorRamp2(c(0, 1, 2 ,3, 4),
                                 c("#0000FFFF", "#9A70FBFF", "#D8C6F3FF", "#FFC8B9FF", "#FF7D5DFF"))
lgd <-  Legend(col_fun = col_fun, 
               title = "log2(count+1)", 
               title_gp = gpar(col="white", cex = 0.75),
               title_position = "leftcenter-rot",
               #direction = "horizontal"
               at = c(0, 1, 4), 
               labels = c("low", "median", "high"),
               labels_gp = gpar(col="white")
               )

绘制图形

1
2
3
4
5

grid.newpage() #新建画布
draw(p) # 绘制热图
draw(lgd, x = unit(0.05, "npc"), 
     y = unit(0.05, "npc"), 
     just = c("left", "bottom")) # 绘制图形

ComplexHeatmap绘制热图非常强大的工具，大部分我想要的功能它都有，甚至我没有想到的它也有，这个教程只是展示其中一小部分功能而已，还有很多功能要慢慢探索。

Tribonacci序列Tn定义:

T0=0, T1=1, T2=1, n>=0时，Tn+3 = Tn + Tn+1 + Tn+2

限制条件是: 0<=n<=37, 32位整型。

我直接用C++撸了下面的代码，

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

#include <iostream>

using namespace std;

class Solution {
public:
    int tribonacci(int n) {
        if (n ==  2){
            return 1;
        } 
        if (n == 1){
            return 1;
        } 
        if (n == 0){
            return 0;
        } 
        int result = tribonacci(n-3) + tribonacci(n-2) + tribonacci(n-1);

        return result;
        
    }
};

// 下面是自己电脑编译测试用的代码
int main(){

	Solution sol;
	for (int i = 0; i <= 37; i++){
	    cout << "tribonacci " << i << " is " << sol.tribonacci(i) << endl;
	}

}

提交答案后，提示”The Limit Exceeded”.

在自己服务器测试时，也发现当n=31之后的速度下降的非常快。

1
2
3

# 编译
g++ -o trib trib.cpp
./trib

原因是当n越大时，直接使用递归会产生大量的重复计算，导致计算效率下降。为了解决该问题，需要维护一个已计算结果的字典，避免重复运算。

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

#include <iostream>

using namespace std;

class Solution {
public:
    int tribo[100];
    int tribonacci(int n) {
        if (n ==  2){
            return 1;
        } 
        if (n == 1){
            return 1;
        } 
        if (n == 0){
            return 0;
        } 
        
        if (tribo[n] > 0){
            return tribo[n]-1;
        }
        
        int result = tribonacci(n-3) + tribonacci(n-2) + tribonacci(n-1);
        tribo[n] = result + 1;

        return result;
        
    }
};

// 下面是自己电脑编译测试用的代码
int main(){

	Solution sol;
	for (int i = 0; i <= 37; i++){
	    cout << "tribonacci " << i << " is " << sol.tribonacci(i) << endl;
	}

}

新的代码中，我预先定义了一个大小为100的数组，主要是因为题目限制n的取值，如果n没有限制取值，我会考虑使用vector容器或者map词典。其次考虑到n=0时，结果为0，并且数组初始化值为0，为了能够通过数组值是否为0来判断是否已经计算了相应的tribonacci值，我对结果加上了一个pseducount。

最后这段代码运行速度是0ms.

如果要用哈希字典，代码如下

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

#include <map>
//一定要引入map头文件
class Solution {
public:
    //定义一个哈希字典
    map<int, int> num;

    int tribonacci(int n) {
        if(n < 2) return n;
        if(n == 2) return 1;
        //如果字典中有已经计算的结果,直接返回结果
        if(num.find(n) != num.end()) return num[n];
        int sum = tribonacci(n-1) + tribonacci(n-2) + tribonacci(n-3);
        num[n] = sum;
        return sum;
    }
};

NextDenovo是武汉未来组胡江博士团队开发的一个三代组装工具，能够用于PacBio和Nanopore数据的组装。但是从官方的介绍而言，此工具在组装Nanopore上优势更大一些。

NextDenovo包括两个模块，NextCorrect用于原始数据纠错，NextGraph能够基于纠错后的进行组装。使用修改版的minimap2进行序列间相互比对。v2.0-beta.1版中在处理高度重复序列上可能存在错误组装，可以通过HiC和BioNano进行纠错。

软件安装

NextDenovo的软件安装非常简单, 下载解压缩即可使用。考虑到NextDenovo需要用Python2.7，我们可以用conda新建一个环境

1
2
3

conda create -n python2 python=2.7
conda activate python2
pip install psutil

然后下载解压缩（我习惯把软件放在~/opt/bisofot下)

1
2
3
4

mkdir -p ~/opt/biosoft/
cd ~/opt/biosoft/
wget https://github.com/Nextomics/NextDenovo/releases/download/v2.0-beta.1/NextDenovo.tgz
tar -zxvf NextDenovo.tgz

测试下软件是否可以使用

1

~/opt/biosoft/NextDenovo/nextDenovo -h

实战

以发表在NC上的拟南芥数据为例, 简单介绍下软件的使用

第一步: 新建一个分析项目

1

mkdir NEXT && cd NEXT

然后从EBI上下载该数据，在run.fofn中记录文件的实际位置。

1
2
3

# 三代测序
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR217/003/ERR2173373/ERR2173373.fastq.gz
realpath ERR2173373.fastq.gz > run.fofn

第二步: 复制和修改配置文件

1

cp ~/opt/biosoft/NextDenovo/doc/run.cfg .

我的配置文件修改如下，参数说明参考官方文档

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

[General]
job_type = local
job_prefix = nextDenovo
task = all # 'all', 'correct', 'assemble'
rewrite = yes # yes/no
deltmp = yes
rerun = 3
parallel_jobs = 5
input_type = raw
input_fofn = input.fofn
workdir = 01_rundir
# cluster_options = -l vf={vf} -q all.q -pe smp {cpu} -S {bash} -w n

[correct_option]
read_cutoff = 1k
seed_cutoff = 3k
blocksize = 3g
pa_correction = 20
seed_cutfiles = 20
sort_options = -m 20g -t 8 -k 40
minimap2_options_raw = -x ava-ont -t 8
correction_options = -p 8

[assemble_option]
random_round = 20
minimap2_options_cns = -x ava-ont -t 8 -k17 -w17
nextgraph_options = -a 1

配置文件的几个重要参数说明（v2.0-beta.1）

• job_type 设置运行环境，可以使用（local， sge， pbs等）
• 运行线程数设置，线程数计算为parallel_jobs分别与sort_option, minimap_options_*的-t数乘积，和correction_options的-p的乘积，量力而行。
• seed_cutfiles 如果在集群上运行，建议设置为可用的节点数，同时设置correction_options的-p为各个节点可用的核数，保证每个节点只有一个correction任务，减少运行时的内存和IO。 如果local上运行， 建议设置为总可用的核除以correction_options的-p值.
• parallel_jobs建议设置至少要大于pa_correction。
• blocksize 是将小于seed_cutfiles的数据拆分成的多个文件时单个文件的大小， 总的比对任务数等于基于该参数切分的文件数乘以seed_cutfiles + seed_cutfiles * (seed_cutfiles - 1)/2， 因此对于10g以内的数据量， 建议设置小于1g， 避免总的任务数小于parallel_jobs的值。
• 测序数据类型相关: 对于PacBio而言，要修改minimap2_options_*中的-x ava-ont为-x ava-pb
• 数据量相关参数: read_cutoff = 1k过滤原始数据中低于1k的read，seed_cutoff = 30k则是选择大于30k以上的数据来矫正。关于seed_cutoff的设置，可以通过~/opt/biosoft/NextDenovo/bin/seq_stat来获取参考值，不建议直接使用默认值，因为改值会受到测序深度和测序长度影响，而且一个不合适的值会显著降低组装质量。对于基因组大于200m以上的物种，-d建议默认。
• correction_options中的-dbuf可以显著降低矫正时的内存，但会显著降低矫正速度。
• random_round参数，建议设置20-100. 该参数是设置随机组装参数的数量，nextGraph会基于每一套随机参数做一次组装， 避免默认参数效果不好。

seq_stat能够根据物种大小和预期用于组装的深度确定seed_cutoff

1

~/opt/biosoft/NextDenovo/bin/seq_stat -g 110Mb -d 30 input.fofn

第三步: 运行NextDenovo

1

~/opt/biosoft/NextDenovo/nextDenovo run.cfg &

运行时间如下

1
2
3

real	64m5.356s
user	1827m37.890s
sys	264m48.246s

默认参数结果是存放在01_rundir/03.ctg_graph/01.ctg_graph.sh.work/ctg_graph00, 可以将其复制到当前目录，用于后续的分析。

1

cat 01_rundir/03.ctg_graph/01.ctg_graph.sh.work/ctg_graph00/nextgraph.assembly.contig.fasta > nextgraph.assembly.contig.fasta

但是在01.ctg_graph.sh.work目录下除了ctg_graph00以外，还有其他随机参数的在组装结果。随机参数结果只输出了统计结果，用户如需要输出组装序列，可以修改01_rundir/03.ctg_graph/01.ctg_graph.sh，将里面的-a 0替换成-a 1。

每个目录下都有shell输出，可以挑选基于nextDenovo.sh.e这里面的结果挑选组装指标较好的，再输出序列，比如说比较下N50

1

grep N50  01_rundir/03.ctg_graph/01.ctg_graph.sh.work/ctg_graph*/*.e

默认情况下，最终组装出20条contig，总大小116M，N50 12M.

使用minimap2将组装结果和比对到TAIR10上，用dotplotly进行可视化

1
2

minimap2 -t 100 -x asm5 Athaliana.fa nextgraph.assembly.contig.fasta > next.paf
dotPlotly/pafCoordsDotPlotly.R  -i next.paf -o next  -l -p 6 -k 5

不难发现，两者存在高度的共线性。大部分TAIR10的染色体对应的都是2条或者3条contig。

此外这篇NC的拟南芥提供了BioNano光学图谱，我使用BioNano Hyrbrid Scaffold 流程进行了混合组装

1
2
3
4
5
6
7
8
9
10

cp /opt/biosoft/Solve3.3_10252018/HybridScaffold/10252018/hybridScaffold_config.xml .
# 修改xml中fasta2cmap的enzyme为BSPQI
perl /opt/biosoft/Solve3.3_10252018/HybridScaffold/10252018/hybridScaffold.pl \
    -n nextgraph.assembly.contig.fasta \
    -b kbs-mac-74_bng_contigs2017.cmap \
    -c hybridScaffold_config.xml \
    -r /opt/biosoft/Solve3.3_10252018/RefAligner/7915.7989rel/RefAligner \
    -o nextgraph \
    -B 2 -N 2 \
    -f

组装结果如下，从原来的20的contig下降到了16个contig。

1
2
3
4
5
6
7

Count  = 16
Min length (Mbp) = 0.026
Median length (Mbp) = 7.224
Mean length (Mbp) = 7.301
N50 length (Mbp) = 13.013
Max length (Mbp) = 14.965
Total length (Mbp) = 116.811

此外还通过BioNano Access进行可视化，以其中一个结果为例。光学图谱和NextDenovo的组装结果存在很高的一致性。

综上，在Nanopore上组装上，我们又多了一个比较好用的工具。

最后非常感谢胡江博士对于本文的指导！

BioNano以cmap格式存放光学图谱，为了评估基因组的组装质量或者了解光学图谱中冗余情况(高杂合基因组组装结果偏大)，我们就需要进行cmap之间的比较。

CMAP间比对

Solve套件提供了runCharacterize.py脚本封装了RefAligner，用于进行CMAP之间的比对。

1
2
3
4
5
6
7

python2.7 runCharacterize.py \
    -t RefAligner的二进制文件路径 \
    -q 用于比对的CMAP \
    -r 参考CMAP \
    -p Pipeline文件路径\
    -a 参数配置文件.xml \
    -n 线程数,默认4

需要注意的是-p和-a参数的设置。-p是Pipeline的文件位置，比如说我的Solve安装在/opt/biosoft/Solve3.4_06042019a，那么参数设置为 -p /opt/biosoft/Solve3.4_06042019a/Pipeline/06042019。 而-a则是要在/opt/biosoft/Solve3.4_06042019a/RefAligner/8949.9232rel/目录下选择合适的xml文件。比如你的CMAP是Irys平台，那么你可以考虑用optArguments_nonhaplotype_irys.xml.

以最新发表的辣椒的光学图谱为例，该物种有比较高的杂合度，组装结果偏大，我们可以通过自比对来寻找冗余区域,

1
2
3
4
5
6
7
8
9

# 下载CMAP
wget https://submit.ncbi.nlm.nih.gov/ft/byid/o62junnn/piper_nigrum_no_rcmap_refinefinal1.cmap
# 自比对
python /opt/biosoft/Solve3.4_06042019a/Pipeline/06042019/runCharacterize.py \
	-t /opt/biosoft/Solve3.4_06042019a/RefAligner/8949.9232rel/RefAligner \
	-q piper_nigrum_no_rcmap_refinefinal1.cmap \
	-r piper_nigrum_no_rcmap_refinefinal1.cmap \
	-p /opt/biosoft/Solve3.4_06042019a/Pipeline/06042019 \
	-a /opt/biosoft/Solve3.4_06042019a/RefAligner/8949.9232rel/optArguments_nonhaplotype_saphyr.xml -n 64

最终会在当前文件下生成一个alignRef文件夹，其中结果是q.cmap,r.cmap和xmap的文件可以用于上传到BioNano Access上进行展示。下图就是一个冗余实例，可以把图中较短的图谱删掉

基因组回帖

为了将基因组回帖到CMAP上，需要先将基因组的fasta格式转成CMAP格式，参数如下

1

perl fa2cmap_multi_color.pl -i 输入FASTA -e 酶1 通道1 [酶2 通道2]

其中一个最重要的参数就是酶切类型。例如我需要将序列回帖到用Nt.BspQI酶切组装的光学图谱上，因此运行参数如下

1

perl /opt/biosoft/Solve3.4_06042019a/HybridScaffold/06042019/scripts/fa2cmap_multi_color.pl -i athaliana.fa -e BspQI 1

最后的athaliana_BSPQI_0kb_0labels.cmap就是模拟酶切的CMAP序列。

之后将模拟酶切的结果回帖到实际的CMAP

1
2
3
4
5
6
7

python /opt/biosoft/Solve3.4_06042019a/Pipeline/06042019/runCharacterize.py \
    -t /opt/biosoft/Solve3.4_06042019a/RefAligner/8949.9232rel/RefAligner \
    -q athaliana_BSPQI_0kb_0labels.cmap \
    -r kbs-mac-74_bng_contigs2017.cmap \
    -p /opt/biosoft/Solve3.4_06042019a/Pipeline/06042019 \
    -a /opt/biosoft/Solve3.4_06042019a/RefAligner/8949.9232rel/optArguments_nonhaplotype_saphyr.xml \
    -n 64

最终会在当前文件下生成一个alignRef文件夹，其中结果是q.cmap,r.cmap和xmap的文件.

SCALE全称是Single-Cell ATAC-seq analysis vie Latent feature Extraction, 从名字中就能知道这个软件是通过隐特征提取的方式分析单细胞ATAC-seq数据。

在文章中，作者从开发者的角度列出了目前的scATAC-seq分析软件，chromVAR, scABC, cisTopic, scVI，发现每个软件都有一定的不足之处，而从我们软件使用者的角度，其实可以考虑都试试这些工具。

SCALE结合了深度生成模型(Depp Generative Models)变分自动编码器框架(Variational Autoencoder, VAE)与概率高斯混合模型(Gaussian Mixture Model, GMM)去学习隐特征，用于准确地鉴定scATAC-seq数据中的特征。

文章通过一张图来解释了软件的工作机制：

SCALE将sc-ATAC-seq的输入数据x(Cells-by-Peaks矩阵)建模成一个联合分布，p(x,z,c)，c是GMM组件中对应的预定义的K个聚类，z是一个隐变量，是细胞在所有peak中实际可能的值，用于后续的聚类和可视化。z通过$z = \mu_z + \sigma_Z \times \epsilon$ 计算而得，公式里面的 $\mu_z$ $\sigma_z$ 是编码器网络从x中学习而得，$\epsilon$ 则是从 $\mathbb{N}(0,\mathbf{I})$ 抽样而成。

从公式中我们还可以发现z其实和GMM的c有关，所以p(x,z,c)也可以写成P(x|z)p(z|c)p(c)，而p(c)是K个预定义聚类分布的离散概率分布，p(z|c)服从混合高斯分布，而p(x|z)则是服从多变量伯努利分布(multivartiable Bernoulli distribution), 通过解码者网络建模而成。

当然从一个软件使用者的角度而言，我们不会去关心代码，也不会关心原理，我们更关心的是这个工具能做什么。SCALE能做以下的分析

• SCALE可以对隐特征聚类识别细胞类群
• SCALE可以降噪，恢复缺失的peak
• SCALE能够区分批次效应和生物学细胞类群之间的差异

软件安装

推荐使用conda的方式进行软件安装（我测试过了，运行没有问题）

第一步：创建一个环境，名字就是SCALE，并且启动该环境

1
2

conda create -n SCALE python=3.6 pytorch
conda activate SCALE

第二步：从GitHub上克隆该项目

1

git clone git://github.com/jsxlei/SCALE.git

第三步：安装SCALE

1
2

cd SCALE
python setup.py install

之后分析的时候，只需要通过conda activate SCALE就能启动分析环境。

考虑后续要交互的读取数据和可视化，那么建议再安装一个Jupyter

1

conda install jupyter

软件使用

SCALE支持两类输入文件：

• count矩阵，行为peak，列为barcode
• 10X输出文件: count.mtx.gz, peak.tsv, barcode.tsv

我们以官方提供的Forebrain数据集为例进行介绍，因为这个数据相对于另外一个数据集Mouse Atlas小多了。

我们在服务器上新建一个文件夹，用于存放从https://cloud.tsinghua.edu.cn/d/bda0332212154163a647/下载的数据

1

mkdir Forebrain

保证Forebrain有下载好的数据

1
2

$ ls Forebrain 
data.txt

之后运行程序

1

SCALE.py -d Forebrain/data.txt -k 8 --impute

软件运行步骤为：

• 加载数据: Loading data
• 模型训练: Training Model
• 输出结果: Saving imputed data

其中模型训练这一步时间比较久，可以尝试用GPU加速（我是普通CPU服务器没有办法）。最终会在当前文件夹看到一个output文件夹，里面有如下内容:

• imputed_data.txt: 每个细胞在每个特征的推断值，建议用--binary保存二进制格式
• model.pt: 用于重复结果的模型文件，--pretrain参数能够读取该模型
• feature.txt: 每个细胞的隐特征，用于聚类和可视化
• cluster_assignments.txt: 两列，barcode和所属类群
• tsne.txt, tsne.pdf: tSNE的坐标和PDF文件，坐标文件可以导入到R语言进行可视化

上面是命令行部分，下面则是Python环境进行交互式操作，输入jupyter notebook，之后在网页上打开

首先是导入各种Python库

1
2
3
4
5
6

import pandas as pd
import numpy as np
from sklearn.metrics import confusion_matrix
from matplotlib import pyplot as plt
import seaborn as sns
from scale.plot import plot_embedding, plot_heatmap

然后加载分析结果，包括聚类信息和特征信息

1
2

y = pd.read_csv('output/cluster_assignments.txt', sep='\t', index_col=0, header=None)[1].values
feature = pd.read_csv('output/feature.txt', sep='\t', index_col=0, header=None)

通过热图展示不同聚类细胞之间的差异图

1
2
3
4
5

plot_heatmap(feature.T, y, 
             figsize=(8, 3), cmap='RdBu_r', vmax=8, vmin=-8, center=0,
             ylabel='Feature dimension', yticklabels=np.arange(10)+1, 
             cax_title='Feature value', legend_font=6, ncol=1,
             bbox_to_anchor=(1.1, 1.1), position=(0.92, 0.15, .08, .04))

如果要矫正批次效应，可以通过根据feature的heatmap，去掉和batch相关的feature来实现

我们可以展示SCALE对原始数据纠正后的值(imputed data), 该结果也能提高chromVAR鉴定motif的效果

1

imputed = pd.read_csv('output/imputed_data.txt', sep='\t', index_col=0)

展示聚类特异性的peak， 分析由mat_specificity_score和cluster_specific完成

1
2
3
4
5
6
7
8

from scale.specifity import cluster_specific, mat_specificity_score

score_mat = mat_specificity_score(imputed, y)
peak_index, peak_labels = cluster_specific(score_mat, np.unique(y), top=200)

plot_heatmap(imputed.iloc[peak_index], y=y, row_labels=peak_labels, ncol=3, cmap='Reds', 
             vmax=1, row_cluster=False, legend_font=6, cax_title='Peak Value',
             figsize=(8, 10), bbox_to_anchor=(0.4, 1.2), position=(0.8, 0.76, 0.1, 0.015))

参数介绍

通过SCALE.py -h可以输出SCALE的所有可用参数

• -d/--dataset: 单个文件矩阵应该指定文件路径，10X输出的多个文件则是文件目录
• -k: 设定输出结果的聚类数
• -o: 输出文件路径
• --pretrain: 读取之前训练的模型
• --lr: 修改起始学习速率, 默认是0.002，和模型训练有关
• --batch_size: 批处理大小， 默认就行，不需要修改（和批次效应处理无关）
• -g GPU: 选择GPU设备数目，非GPU服务器用不到
• --seed: 初始随机数种子，通常在遇到nan缺失时考虑修改
• -encode_dim, -decode_dim: 编码器和解码器的维度，通常也不需要修改
• -latent 隐藏层维度
• --low, --high: 过滤低质量的peak, 即出现比例高于或者低于某个阈值的peak，默认是0.01和0.9。作者推荐保留1万-3万的peak用于SCALE分析。如果数据质量很高，且peak数不多于10万，那么可以不过滤。
• --min_peaks: 过滤低质量细胞，如果该细胞的peak低于阈值
• log_transform: log2(x+1)的变换
• --max_iter: 最大迭代数，默认是30000, 可以观察损失收敛的情况来修改，也就是训练模型这一步输出的信息
• -weight_decay: 没有说明
• --impute: 保存推断数据，默认开启
• --binary: 推荐加上该参数，减少imputed data占用空间
• --no_tsne: 不需要保存t-SNE结果
• --reference: 参考细胞类型
• -t: 如果输出矩阵是列为peak，行为barcode，用该参数进行转置

对于使用者而言，我们一般只用修改-k更改最后的聚类数，--low, --high, ---min_peaks来对原始数据进行过滤，以及加上--binary节约空间。

假如在训练模型阶段，发现输出信息为loss=nan recon_loss=nan kl_loss=nan,十有八九最终会报错退出， 可以如下的参数调整

• 更改--seed
• 用更加严格的条件过滤peak，例如-x 0.04 或 -x 0.06
• 降低初始的学习速率，--lr 0.0002

文献出处: Sequencing a Juglans regia × J. microcarpa hybrid yields high-quality genome assemblies of parental specie

文章的亮点在于通过对一个F1子代进行三代测序，之后利用BioNano组装出两个亲本的光学图谱，最后根据光学图谱从F1中拆分出两套单倍型。

杂合基因组组装的时候，通常会通过构建近交系，花粉加倍(DH)来降低杂合度。文章采用了种间杂交的方式来避免杂合。这是因为种间杂交的基因组通常包括亲本的两个单倍型，因此能够从中分离出两个亲本的单倍型基因组。种间杂交的一个优势在于，它比较容易构建。

分析策略流程图如下

文字版讲解如下：

Step1: 使用两个近源物种（都有基因组草图）进行杂交得到F1(MS1-56)，使用PacBio进行测序组装（460条contig, 1.05G, 两倍基因组大小）

Step2: 使用BioNano构建F1和两个亲本的光学图谱。

Step3: 通过自我比对（self-alignment），鉴定haploid和diploid。 亲本的光学图谱包含两个部分，一个是haploid（两个单倍体因为相似而坍缩成一个)，另一个是diploid(两个单倍体因为不够相似被拆分， 也就是分型(phased))。大概是下面这个效果，a表示Serr ctg4上存在两套单倍型的情况.

Step4: 将contig回帖到杂合的光学图谱上，组装成scaffolds，大小为1.06G。之后根据杂合后代光学图谱和亲本光学图谱的关系，对scaffold进行拆分。上图b,c表示的是将F1后代的光学图谱和双亲的光学图谱进行比较，来鉴定出子代中属于一方亲本的单倍型。

Step5: 上述因为太短而无法回帖到光学图谱的contig，比对到illumina组装的双亲基因组，进行区分。

Step6: 根据高密度遗传图谱对scaffold进行排序和确认方向，最终得到准染色体级别的基因组。

最后基因组还可以补洞和纠错，来提高质量，比如说10X Genomics Chicago 测序。

最后一点感想：文章利用光学图谱进行分型的思路的确很棒。但是我想的是，为啥不直接测亲本基因组，用HiC进行分型呢？

C++能解决的瓶颈问题有:

• 由于迭代依赖于之前结果，循环难以简便的向量化运算
• 递归函数，或者是需要对同一个函数运算成千上万次
• R语言缺少一些高级数据结构和算法

我们只需要在代码中写一部分C++代码来就可以处理上面这些问题。后续操作在Windows下进行，你需要安装Rtools，用install.packages("Rcpp")安装新版的Rcpp，最重要一点，你需要保证你R语言时不能是C:/Program Files/R/R-3.5.1/这种形式，否则会报错。

后续操作会用到microbenchmark包来评估R代码和RCPP的效率差异，用install.packages('microbenchmark)安装

RCPP入门

先从一个简单的add函数开始，学习如何用cppFunction在R里面写C++代码

1
2
3
4
5
6
7
8
9

library(Rcpp)

cppFunction('int add(int x, int y, int z) {
  int sum = x + y + z;
  return sum;
}')
add
# function (x, y) 
# .Call(<pointer: 0x0000000063c015a0>, x, y)

Rcpp将会编译C++代码， 然后构建能够连接到C++函数的R函数。后续将会介绍如何将一些R代码改写成C++代码。

• 标量输入，标量输出
• 向量输入，标量输出
• 向量输入，向量输出
• 矩阵输入，向量输出

没有输入，标量输出

最简单的函数就是不提供任何输出，返回一个输出，比如说

1

one <- function() 1L

等价的C++代码是

1
2
3

int one(){
    return 1;
}

那么将这段C++代码在R用cppFunction中改写就是如下

1
2
3

cppFunction('int one(){
  return 1;
}')

上面这段函数就展示了R和C++之间一些重要区别:

• C++写代码不是 函数名 <- function(参数){} 而是 函数名(函数参数){}
• C++中必须声明返回类型，ini就是标量整数。C++对应R语言常用向量的类是: NumericVector,IntegerVector, CharacterVector和LogicalVector.
• R语言没有标量，全是向量。而C++有向量和标量之分，标量的数据类型是double, int, String和bool
• C++你必须要用到return声明要返回的数据
• 每段代码后要跟着;

标量输入，标量输出

我们可以写一个函数，sign，他的功能就是把一个负数转成正数，正数不变

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

signR <- function(x){
  if (x > 0){
    x
  } else if (x == 0 ){
    0
  } else{
    -x
  }
}

cppFunction('int signC(int x){
            if( x >0 ){
              return x;
            } else if (x == 0){
              return 0;
            } else {
              return -x;
            }
}')

这个例子中要注意两件事情

• 在C++中，你需要声明输入的数据类型
• C++和R的条件语句长得一样。

向量输入，标量输出

R和C++一大区别就是R的循环效率很低。因此在R语言要尽量避免使用显示的循环语句，尽量向量化运算函数。而C++的循环花销特别小，所以可以放心大胆的用。

让我们用R代码写一个求和函数sum 以及 C++的求和函数，然后比较下效率

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

sumR <- function(x){
  total <- 0
  for (i in seq_along(x)){
    total <- total + x[i]
  }
  total
}

cppFunction('int sumC(NumericVector x ){ 
            int n = x.size();
            double total = 0;
            for(int i = 0; i < n; ++i){
              total += x[i];  
            }
            return total;
            }')

C++版本和R版本的逻辑相同，但是有如下不同

• 用.size()确认向量的长度
• for的写法为for(初始值; 判断语句; 递增)
• 记住: C++的向量索引从0开始，R是从1开始
• 向量赋值是=而不是<-
• total += x[i]等价于total = total + x[i], 类似的符号还有-=, *=, /=

最后用microbenchmark比较下,R自带求和函数和我们自己写的两个版本的差异

1
2
3
4
5
6

x <- runif(1000)
microbenchmark(
  sum(x),
  sumC(x),
  sumR(x)
)

最快的是高度优化过的内置函数，最差的就是sumR(), 速度会比sumC()慢10倍以上。

向量输入，向量输出

R中比较常见的操作就是向量间运算，尤其R还会自动补齐。自动补齐某些时候会造成一些问题，但是C++不存在这个问题。我们可以写一个RCPP的+函数

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

cppFunction('NumericVector addC(NumericVector x, NumericVector y){
  int xn = x.size();
  int yn = y.size();
  
  if (xn != yn){
    stop("input should be same length");
  }
  NumericVector out(xn);
  for(int i=0; i< xn; ++i){
    out[i] = x[i] + y[i];
  }
  return out;
}')

x <- runif(1e6)
y <- runif(1e6)
microbenchmark(addC(x,y),
               x+y)

矩阵输入，向量输出

每个向量类型都有矩阵等价类，NumericMatrix, IntegerMatirx, CharacterMatirx, LogicalMatirx. 让我们尝试写一个rowSums()函数

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

cppFunction('NumericVector rowSumsC(NumericMatrix x){
  int nrow = x.nrow(), ncol = x.ncol();
  NumericVector out(nrow);
  
  for(int i = 0; i < nrow; i++){
    double total =0;
    for(int j =0; j< ncol; j++){
      total += x(i,j);
    }
    out[i] = total;
  }
  return out;
}')
set.seed(1024)
x <- matrix(sample(100), nrow = 10)
rowSumsC(x)

这里注意有两点不同，在C++中，你用()对矩阵取值，而不是[]

尽管看起来C++的代码运行起来比R语言快多了，比如说R要一分钟，RCPP只要一秒，但是如果算上我们写代码的时间和调试代码的时间，刚开始不熟练估计要10分钟，那么总体来看，还是直接上手写R代码比较合适。

但是如果有一些代码要不断复用，那么写C++代码还是很划算。这个时候就建议将代码写到专门的文本中，用sourceCpp()加载，而不是cppFunction()函数

在Rsutdio中可以创建一个C++模板文件，代码写完之后还可以进行debug。

比如说在里面写上面的rowSumsC函数，分为如下几个部分

导入头文件，加载Rcpp到命名空间中，类似于library()

1
2

#include <Rcpp.h>
using namespace Rcpp;

使用// [[Rcpp::export]]说明这里的函数会被R使用

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

// [[Rcpp::export]]
NumericVector rowSumsC(NumericMatrix x){
  int ncol = x.ncol(), nrow = x.nrow();
  NumericVector out(nrow);
  
  for (int i =0; i < nrow; i++ ){
    double total = 0;
    for (int j =0 ;j < ncol; j++){
      total += x(i,j);
    }
    out[i] = total;
  }
  return out;
}

下面部分会在sourceCpp()加载后自动运行

1
2
3
4
5
6
7
8
9

/*** R
library(microbenchmark)
set.seed(1014)
x <- matrix(sample(100), 10)
microbenchmark(
  rowSumsC(x),
  Matrix::rowSums(x)
)
*/

将文件保存成rowSumsC.cpp, 之后在R里用sourceCpp(file = "rowSumsC.cpp")。