使用Swiss-Prot根据同源基因进行注释

第一步: 在uniprot下载UniProt 上植物dat格式的注释文件。

1
2
wget ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/taxonomic_divisions/uniprot_sprot_plants.dat.gz
wget ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/taxonomic_divisions/uniprot_trembl_plants.dat.gz

将两个dat合并到成一个文件

1
zcat uniprot_sprot_plants.dat.gz uniprot_trembl_plants.dat.gz > uniprot_plants.dat

第二步: 从dat中提取fasta序列

1
2
dat=uniprot_plants.dat
awk '{if (/^ /) {gsub(/ /, ""); print} else if (/^AC/) print ">" $2}' $1 > ${1%%.dat}_AC.fasta

第三步: 建立DIAMOND或NCBI BLAST+索引

1
diamond makedb --in uniprot_plants_AC.fasta -d uniprot_plants_AC

第四步: 使用DIAMOND或NCBI BLAST+进行比对

1
diamond blastp -d /data/database/UniProt-Plant/uniprot_plants_AC.dmnd -q proteins.fasta --evalue 1e-5 > blastp.outfmt6

第五步: 从DIMAMOND或NCBI BLAST+的比对结果中筛选每个query的最佳subject

1
python -m jcvi.formats.blast best -n 1 blastp.outfmt6

第六步: 使用add_annotation_from_dat.py(代码在GitHub上)根据blastp输出从dat中提取GO/KEGG/同源基因。运行在Python2/3环境中,需要安装BioPython

1
python ~/myscripts/add_annotation_from_dat.py blastp.outfmt6.best /data/database/UniProt-Plant/uniprot_plants.dat

之后会输出swiss_annotation.tsv, 输出信息包括如下几列

  • gene id
  • uniprot accession
  • identity
  • homology species
  • EnsemblPlants
  • GO ID
  • GO component, CC/MF/BP
  • evidence

NECAT:Nanopore数据的高效组装工具

NECAT是肖传乐老师团队开发的一个针对Nanopore数据组装的软件,目前该工具尚未发表,除了https://github.com/xiaochuanle/NECAT有软件的介绍外,暂时没有中文资料介绍NECAT的使用。

太长不看的结论: Nanopore的组装推荐用下NECAT。组装之后是先用MEDAKA做一遍三代polish,然后用NextPolish默认参数做二代polish。

这篇将会以一篇发表在Nature Communication上的拟南芥nanopore数据介绍如何使用NECAT进行组装,运行在CentOS Linux release 7.3.1611 (Core),64G为内存, 20线程(Intel(R) Xeon(R) CPU E5-2640 v4 @ 2.40GHz),下面是正文。

软件安装

NECAT可以在https://github.com/xiaochuanle/NECAT/releases/页面获取最新的软件下载地址,这里下载的是0.01版本。

1
2
3
wget https://github.com/xiaochuanle/NECAT/releases/download/v0.01/necat_20190307_linux_amd64.tar.gz
tar xzvf necat_20190307_linux_amd64.tar.gz
export PATH=$PATH:$(pwd)/NECAT/Linux-amd64/bin

目前0.01版本不支持gz文件作为输入,但后续版本应该会支持。

目前更新到necat_20200119_Linux-amd64,新版本安装方法为

1
2
3
4
wget https://github.com/xiaochuanle/NECAT/releases/download/SourceCodes20200119/necat_20200119_Linux-amd64.tar.gz
tar xzvf necat_20200119_Linux-amd64.tar.gz
cd NECAT/Linux-amd64/bin
$ export PATH=$PATH:$(pwd)

新版本增加gz文件支持。目前测试发现文件名需要符合xxx.fastqxxx.fastq.gz命名格式,对于fq.gz无法识别,会导致程序文件出错。

实战

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

1
mkdir NECAT && cd NECAT

以发表在NC上的拟南芥数据为例, 下载该数据

1
2
3
# 三代测序
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR217/003/ERR2173373/ERR2173373.fastq.gz
seqkit seqkit fq2fa ERR2173373.fastq.gz | gzip -c > ERR2173373.fasta

第二步: 创建配置文件

1
necat.pl config ath_config.txt 

配置文件中,主要修改如下几个参数

1
2
3
4
5
6
PROJECT=athaliana #项目名
ONT_READ_LIST=read_list.txt #read所在路径文件
GENOME_SIZE=120000000 #基因组大小
THREADS=20 # 线程数
MIN_READ_LENGTH=3000 # 最短的read长度
CNS_OUTPUT_COVERAGE=45 # 用于组装的深度

参数中还有一个,NUM_ITER=2,它并非是简单的重复2次纠错,它的每一轮的校正目的其实不同,第一轮的优先级是敏感度(senstitive), 第二轮之后主要追求速度(fast)。

除了上面的配置参数外,其他参数可以不需要修改,使用默认的值即可。需要修改的话,参考最后的参数说明部分。

第三步: 序列纠错

1
necat.pl correct ath_config.txt &

纠错后的reads在athaliana/1-consensus/cns_final.fasta

cns_finla.fasta的统计信息会输出在屏幕中, 或者自己用fsa_rd_stat也能得到同样的结果

1
2
3
4
5
6
7
8
9
10
Count: 206342
Tatal: 3102480870
Max: 112992
Min: 1010
N25: 31940
L25: 18989
N50: 21879
L50: 48506
N75: 13444
L75: 93215

此外我还用time获取了运行时间,纠错花了大概一个小时。

1
2
3
real	55m31.451s
user 815m32.801s
sys 7m55.039s

第四步: contig组装

1
necat.pl assemble ath_config.txt &

结果在athaliana/4-fsa/contigs.fasta

关于contigs.fata统计信息会输出在屏幕上,同样用fsa_rd_stat 也可以。

1
2
3
4
5
6
7
8
9
10
Count: 162
Tatal: 122293198
Max: 14562810
Min: 1214
N25: 13052494
L25: 3
N50: 9503368
L50: 5
N75: 4919866
L75: 10

时间用了75分钟

1
2
3
real	74m53.127s
user 1308m29.534s
sys 12m5.032s

第五步: contig搭桥

1
necat.pl bridge ath_config.txt

结果在athaliana/6-bridge_contigs/bridged_contigs.fasta

1
2
3
4
5
6
7
8
9
10
Count: 127
Tatal: 121978724
Max: 14562810
Min: 2217
N25: 13193939
L25: 3
N50: 11146374
L50: 5
N75: 5690371
L75: 9

从N50和N75可以看出这一步会提高组装的连续性。

组装结果polish

对Nanopore组装结果进行polish的常用软件有下面3个

由于拟南芥的基因组比较小,我分别用了Medaka和racon对输出结果进行polish(因为没有原始信号数据,因此nanopolish用不了),代码如下

Medaka

1
2
3
4
5
NPROC=20
BASECALLS=ERR2173373.fasta
DRAFT=athaliana/6-bridge_contigs/bridged_contigs.fasta
OUTDIR=medaka_consensus
medaka_consensus -i ${BASECALLS} -d ${DRAFT} -o ${OUTDIR} -t ${NPROC} -m r941_min_high

三轮Racon:

1
2
3
4
5
6
7
gzip -dc ERR2173373.fastq.gz  > ERR2173373.fastq
minimap2 -t 20 ${DRAFT} ERR2173373.fastq > round_1.paf
racon -t 20 ERR2173373.fastq round_1.paf ${DRAFT} > racon_round1.fasta
minimap2 -t 20 racon_round1.fasta ERR2173373.fastq > round_2.paf
racon -t 20 ERR2173373.fastq round_2.paf racon_round1.fasta> racon_round2.fasta
minimap2 -t 20 racon_round2.fasta ERR2173373.fastq > round_3.paf
racon -t 20 ERR2173373.fastq round_3.paf racon_round2.fasta> racon_round3.fasta

在后续评估质量的时候,我发现单纯用三代polish的结果还不是很好,因此我用他们提供的二代测序,用NextPolish对NECAT的结果进行polish。

1
2
3
# 二代测序
prefetch ERR2173372
fasterq-dump -O . ERR2173372

run.cfg内容如下, 其中sgs.fofn记录的就是解压后的ERR2173372_1.fastq和ERR2173372_2.fastq的路径

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
[General]                
job_type = local
job_prefix = nextPolish
task = 1212
rewrite = no
rerun = 3
parallel_jobs = 8
multithread_jobs = 20
genome = input.fasta
genome_size = auto
workdir = ./nextpolish
polish_options = -p {multithread_jobs}
[sgs_option]
sgs_fofn = ./sgs.fofn
sgs_options = -max_depth 100 -bwa

我考虑了两种情况,一种是直接用二代polish,另一种是三代polish之后接二代polish。

结果评估

计算时间上,我之前用Canu跑了相同的数据,设置原始错误率0.5,纠错后错误率为0.144,用3个节点(每个节点12个线程),运行了3天时间,但是NECAT只需要3个小时左右就能完成相同的分析,这个速度差异实在是太明显了。

用Minimap2 + dotPlotly绘制CANU,NECAT和拟南芥参考基因组的共线性图

1
2
3
4
minimap2 -t 20 -x asm5 Athaliana.fa NECAT.fa > NECAT.paf
pafCoordsDotPlotly.R -i NECAT.paf -o NECAT -l -p 10 -k 5
minimap2 -t 20 -x asm5 Athaliana.fa CANU.fa > CANU.paf
pafCoordsDotPlotly.R -i CANU.paf -o CANU -l -p 10 -k 5

NECAT的结果

NECAT

CANU的结果

CANU

NECAT和CANU都和参考基因组有着良好的共线性,但是NECAT的连续性更好,几乎成一条直线。

之后,我使用了QUAST来评估Canu,NECAT初步组装,NECAT用Medaka, nanopolish和racon纠错的结果(MD: MEDAKA, RC: RACON, NP:NextPolish)。

1
2
3
4
5
6
7
8
9
10
quast.py -t 100 --output-dir athaliana --circos \
CANU.fa \
NECAT.fa \
NECAT_MD.fa \
NECAT_MD_NP.fa \
NECAT_NP.fa \
NECAT_RC.fa \
NECAT_RC_NP.fa \
-r Athaliana.fa \
-g TAIR10_GFF3_genes.gff &

一些描述基本信息

1
2
3
4
5
6
7
CANU        N50 = 4875070,  L50 = 7, Total length = 114689024, GC % = 36.09 
NECAT N50 = 11146374, L50 = 5, Total length = 121978724, GC % = 36.50
NECAT_MD N50 = 11216803, L50 = 5, Total length = 122101599, GC % = 36.54
NECAT_MD_NP N50 = 11405151, L50 = 5, Total length = 124142955, GC % = 36.30
NECAT_NP N50 = 11399084, L50 = 5, Total length = 124735066, GC % = 36.36
NECAT_RC N50 = 11212098, L50 = 5, Total length = 122519370, GC % = 36.4
NECAT_RC_NP N50 = 11406553, L50 = 5, Total length = 124618502, GC % = 36.34

在BUSCO完整度上, 以embryophyta_odb10作为物种数据库, 其中ONTmin_IT4是发表的文章里的结果, Athalina则是拟南芥的参考基因组,我们以它们的BUSCO值作为参照。

1
2
3
4
5
6
Athalina     : C:98.6%[S:98.0%,D:0.6%],F:0.4%, M:1.0%, n:1375
ONTmin_IT4 : C:98.4%[S:97.7%,D:0.7%],F:0.7%, M:0.9%, n:1375
CANU : C:22.9%[S:22.8%,D:0.1%],F:20.2%,M:56.9%,n:1375
NECAT : C:36.6%[S:36.6%,D:0.0%],F:22.9%,M:40.5%,n:1375
NECAT_MEDAKA : C:53.6%[S:53.2%,D:0.4%],F:21.0%,M:25.4%,n:1375
NECAT_RACON : C:45.3%[S:45.2%,D:0.1%],F:23.1%,M:31.6%,n:1375

二代Polish后的BUSCO结果如下(MD: MEDAKA, RC: RACON, NP:NextPolish):

1
2
3
4
5
Athalina   : C:98.6%[S:98.0%,D:0.6%],F:0.4%,M:1.0%,n:1375
ONTmin_IT4 : C:98.4%[S:97.7%,D:0.7%],F:0.7%,M:0.9%,n:1375
NECAT_NP : C:98.6%[S:97.9%,D:0.7%],F:0.4%,M:1.0%,n:1375
NECAT_MD_NP: C:98.7%[S:98.0%,D:0.7%],F:0.4%,M:0.9%,n:1375
NECAT_RC_NP: C:98.5%[S:97.8%,D:0.7%],F:0.4%,M:1.1%,n:1375

从以上这些数据,你可以得到以下几个洞见:

  • 在Nanopore的组装上,NECAT效果优于Canu,无论是连续性还是N50上
  • MEDAKA三代polish效果好于RACON。在速度上,MEDAKA比三遍RACON都慢,并且MEDAKA会将一些可能的错误组装给打断
  • Nanopore的数据用NECAT组装后似乎用NextPolish进行polish后就行,但是由于物种比较小,可能不具有代表性。

结论: Nanopore的组装建议用NECAT。组装之后是先用MEDAKA做一遍三代polish,然后用NextPolish默认参数做二代polish。

配置文件补充

这部分对配置文件做一点简单补充。

下面这些参数相对简单,不需要过多解释,按照自己需求修改

  • CLEANUP: 运行完是否清理临时文件,默认是0,表示不清理
  • USE_GRID: 是否使用多节点, 默认是false
  • GRID_NODE: 使用多少个节点,默认是0,当USE_GRID为true时,按照自己实际情况设置

以下的参数则是需要根据到具体的软件中去查看具体含义,需要和软件开发者讨论

  • OVLP_FAST_OPTIONS: 第二轮纠错时, 传给oc2pmov
  • OVLP_SENSITIVE_OPTIONS: 第一轮纠错时, 传给oc2pmov
  • CNS_FAST_OPTIONS: 第二轮纠错时,传给oc2cns
  • CNS_SENSITIVE_OPTIONS: 第一轮纠错时,传给oc2cns
  • TRIM_OVLP_OPTIONS: 传给oc2asmpm
  • ASM_OVLP_OPTIONS: 传给oc2asmpm
  • FSA_OL_FILTER_OPTIONS: 参数传给fsa_ol_filter
  • FSA_ASSEMBLE_OPTIONS: 参数传给fsa_assemble
  • FSA_CTG_BRIDGE_OPTIONS: 参数传给fsa_ctg_bridge

参考资料

Seurat的scATAC-seq分析流程

Seurat 3.X版本能够整合scRNA-seq和scATAC-seq, 主要体现在:

  • 基于scRNA-seq的聚类结果对scATAC-seq的细胞进行聚类
  • scRNA-seq和scATAC-seq共嵌入(co-embed)分析

整合步骤包括如下步骤:

  1. 从ATAC-seq中估计RNA-seq表达水平,即从ATAC-seq reads定量基因表达活跃度
  2. 使用LSI学习ATAC-seq数据的内部结构
  3. 鉴定ATAC-seq和RNA-seq数据集的锚点
  4. 数据集间进行转移,包括聚类的标签,在ATAC-seq数据中推测RNA水平用于共嵌入分析

数据下载

测试数据下载地址:

可以复制下载链接到浏览器下载,也可以直接在R语言用download.file中进行下载。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# 下载peak
atac_peak <- "http://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5"
atac_peak_file <- "atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5"
download.file(atac_peak, atac_peak_file)

# 下载singlecell.csv
singlecell <- "http://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_singlecell.csv"
singlecell_file <- "atac_v1_pbmc_10k_singlecell.csv"
download.file(singlecell, singlecell_file)

# 下载rna-seq
rna_seq <- "http://cf.10xgenomics.com/samples/cell-exp/3.0.0/pbmc_10k_v3/pbmc_10k_v3_filtered_feature_bc_matrix.h5"
rna_seq_file <- "pbmc_10k_v3_filtered_feature_bc_matrix.h5"
download.file(rna_seq, rna_seq_file)

# 下载GTF
gtf <- "ftp://ftp.ensembl.org/pub/grch37/release-82/gtf/homo_sapiens/Homo_sapiens.GRCh37.82.gtf.gz"
gtf_file <- "Homo_sapiens.GRCh37.82.gtf.gz"
download.file(gtf, gtf_file)

基因活跃度定量

首先,先将peak矩阵转成基因活跃度矩阵。Seurat做了一个简单的假设,基因活跃度可以通过简单的将落在基因区和其上游2kb的count相加得到基因活跃度,并且这个结果Cicero等工具返回gene-by-cell矩阵是类似的。

1
2
3
4
5
6
7
8
9
library(Seurat)
library(ggplot2)
# 读取peak
peaks <- Read10X_h5(atac_peak_file)
activity.matrix <- CreateGeneActivityMatrix(peak.matrix = peaks,
annotation.file = gtf_file,
seq.levels = c(1:22, "X", "Y"),
upstream = 2000,
verbose = TRUE)

activity.matrix是一个dgCMatrix对象,其中行为基因,列为细胞。因此如果对于Cicero的输出结果,只要提供相应的矩阵数据结构即可。

设置对象

下一步,我们要来设置Seurat对象,将原始的peak counts保存到assay中,命名为”ATAC”

1
pbmc.atac <- CreateSeuratObject(counts = peaks, assay = "ATAC", project = "10x_ATAC")

增加基因活跃度矩阵,命名为”ACTIVITY”.

1
pbmc.atac[["ACTIVITY"]] <- CreateAssayObject(counts = activity.matrix)

增加细胞的元信息,该信息来自于scATAC-seq的CellRanger处理的singlecell.csv

1
2
3
4
meta <- read.table(singlecell_file, sep = ",", header = TRUE, row.names = 1, 
stringsAsFactors = FALSE)
meta <- meta[colnames(pbmc.atac), ]
pbmc.atac <- AddMetaData(pbmc.atac, metadata = meta)

过滤掉scATAC-seq数据中总count数低于5000的细胞,这个阈值需要根据具体实验设置

1
2
pbmc.atac <- subset(pbmc.atac, subset = nCount_ATAC > 5000)
pbmc.atac$tech <- "atac"

数据预处理

为了找到scATAC-seq数据集和scRNA-seq数据集之间的锚定点,我们需要对基因活跃度矩阵进行预处理

设置pbmc.atac的默认Assay为”ACTIVITY”, 然后寻找高表达的基因,对基因活跃度矩阵进行标准化和Scale。

1
2
3
4
DefaultAssay(pbmc.atac) <- "ACTIVITY"
pbmc.atac <- FindVariableFeatures(pbmc.atac)
pbmc.atac <- NormalizeData(pbmc.atac)
pbmc.atac <- ScaleData(pbmc.atac)

同样的,我们还需要对peak矩阵进行处理。这里我们用的隐语义(Latent semantic indexing, LSI)方法对scATAC-seq数据进行降维。该步骤学习scRNA-seq数据的内部结构,并且在转换信息时对锚点恰当权重的决定很重要。

根据 Cusanovich et al, Science 2015提出的LSI方法,他们搞了一个RunLSI函数。LSI的计算方法为TF-IDF加SVD。

我们使用在所有细胞中至少有100个read的peak,然后降维到50。该参数的选择受之前的scATAC-seq研究启发,所以没有更改,当然你你可以把它改了。

1
2
3
4
5
DefaultAssay(pbmc.atac) <- "ATAC"
VariableFeatures(pbmc.atac) <- names(which(Matrix::rowSums(pbmc.atac) > 100))
pbmc.atac <- RunLSI(pbmc.atac, n = 50, scale.max = NULL)
#pbmc.atac <- RunUMAP(pbmc.atac, reduction = "lsi", dims = 1:50)
pbmc.atac <- RunTSNE(pbmc.atac, reduction = "lsi", dims = 1:50)

: 要将pbmc.atac的默认assay切换成”ATAC”, 非线性降维可以选择UMAP或者t-SNE。

我们之前使用过Seurat对scRNA-seq数据进行预处理和聚类,下载地址为Dropbox

1
2
pbmc.rna <- readRDS("pbmc_10k_v3.rds")
pbmc.rna$tech <- "rna"

将scRNA-seq和scATAC-seq共同展示,对一些骨髓(myeloid)和淋巴(lymphoid)PBMC中比较常见的聚类,其实是能从图中看出来。

1
2
3
4
5
6
7
p1 <- DimPlot(pbmc.atac, reduction = "tsne") + 
NoLegend() +
ggtitle("scATAC-seq")
p2 <- DimPlot(pbmc.rna, group.by = "celltype", reduction = "tsne", label = TRUE, repel = TRUE) +
NoLegend() +
ggtitle("scRNA-seq")
CombinePlots(plots = list(p1, p2))

标签转移前

现在,我们用FindTransferAnchors鉴定scATAC-seq数据集和scRNA-seq数据集的锚点,然后根据这些锚点将 10K scRNA-seq数据中鉴定的细胞类型标记转移到scATAC-seq细胞中。

1
2
3
4
5
6
transfer.anchors <- FindTransferAnchors(reference = pbmc.rna, 
query = pbmc.atac,
features = VariableFeatures(object = pbmc.rna),
reference.assay = "RNA",
query.assay = "ACTIVITY",
reduction = "cca")

Seurat比较的是scRNA-seq表达量矩阵和scATAC-seq中基因活跃度矩阵,利用CCA降维方法比较两者在scRNA-seq中的高变异基因的关系

为了转移细胞类群的编号,我们需要一组之前注释过的细胞类型,作为TransferData的 refdata 参数输入。TransferData本质上采用的是KNN算法,利用距离未知类型细胞的最近N个已知细胞所属的类型来定义该细胞。weight.reduction参数是用来选择设置权重的降维。

1
2
3
4
celltype.predictions <- TransferData(anchorset = transfer.anchors, 
refdata = pbmc.rna$celltype,
weight.reduction = pbmc.atac[["lsi"]])
pbmc.atac <- AddMetaData(pbmc.atac, metadata = celltype.predictions)

我们可以检查每个细胞预测得分的分布情况,选择性的过滤哪些得分比较低的细胞。我们发现超过95%的细胞得分大于或等于0.5.

1
2
hist(pbmc.atac$prediction.score.max)
abline(v = 0.5, col = "red")

得分分布

1
2
3
table(pbmc.atac$prediction.score.max > 0.5)
# FALSE TRUE
# 383 7483

之后,我们就可以在UMAP上检查预测的细胞类型的分布,检查细胞类型在scRNA-seq和scATAC-seq中的相对位置

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
pbmc.atac.filtered <- subset(pbmc.atac, 
subset = prediction.score.max > 0.5)
pbmc.atac.filtered$predicted.id <- factor(pbmc.atac.filtered$predicted.id,
levels = levels(pbmc.rna)) # to make the colors match
p1 <- DimPlot(pbmc.atac.filtered,
group.by = "predicted.id",
label = TRUE,
repel = TRUE) +
ggtitle("scATAC-seq cells") +
NoLegend() +
scale_colour_hue(drop = FALSE)
p2 <- DimPlot(pbmc.rna, group.by = "celltype", label = TRUE, repel = TRUE) +
ggtitle("scRNA-seq cells") +
NoLegend()
CombinePlots(plots = list(p1, p2))

预测后结果

在转移细胞类型标记之后,我们就可以进行细胞特意水平上的下有分析。举个例子,我们可以去找一些某些细胞类型特异的增强子,寻找富集的motif。目前这些分析Seurat还不直接支持,还在调试中。

共嵌入(co-embedding)

最后,如果你想将所有的细胞一同展示,可以将scRNA-seq和scATAC-seq数据嵌入到相同的低维空间。

我们使用之前的锚点从scATAC-seq细胞中推断RNA-seq的值,后续分析就相当于两个单细胞数据的分析流程。

注意: 这一步只是为了可视化,其实不做也行。

选择用于估计的基因,可以是高变动基因,也可以是所有基因。

1
2
3
# 只选择高变动的基因作为参考
genes.use <- VariableFeatures(pbmc.rna)
refdata <- GetAssayData(pbmc.rna, assay = "RNA", slot = "data")[genes.use, ]

之后利用TransferData推断scATAC-seq在这些基因中的可能值,输出结果就是ATAC细胞对应的scRNA-seq矩阵

1
2
3
4
5
imputation <- TransferData(anchorset = transfer.anchors, 
refdata = refdata,
weight.reduction = pbmc.atac[["lsi"]])
# this line adds the imputed data matrix to the pbmc.atac object
pbmc.atac[["RNA"]] <- imputation

合并两个的结果,然后就是scRNA-seq的常规分析。

1
2
3
4
5
6
7
coembed <- merge(x = pbmc.rna, y = pbmc.atac)
# 标准化分析
coembed <- ScaleData(coembed, features = genes.use, do.scale = FALSE)
coembed <- RunPCA(coembed, features = genes.use, verbose = FALSE)
#coembed <- RunUMAP(coembed, dims = 1:30)
coembed <- RunTSNE(coembed, dims = 1:30)
coembed$celltype <- ifelse(!is.na(coembed$celltype), coembed$celltype, coembed$predicted.id)

在t-SNE上绘制结果

1
2
3
p1 <- DimPlot(coembed, reduction="tsne", group.by = "tech")
p2 <- DimPlot(coembed, reduction="tsne", group.by = "celltype", label = TRUE, repel = TRUE)
CombinePlots(list(p1, p2))

tSNE plot

从上面的结果中,你可能会发现某些细胞只有在一类技术中存在。举个例子,从巨噬细胞(megakaryocytes)到成熟的血小板细胞(pletelet)由于没有细胞核,无法被scATAC-seq技术检测到。我们可以单独展示每个技术,进行检查

1
DimPlot(coembed, reduction="tsne", split.by = "tech", group.by = "celltype", label = TRUE, repel = TRUE) + NoLegend()

分别展示

此外,你还可以发现B细胞前体类型附近有一类细胞只由scATAC-seq的细胞构成。通过展示这些细胞在CellRanger分析结果提供的黑名单位置所占read数,可以推测出这类细胞可能是死细胞,或者是其他scRNA-seq无法检测的人为因素。

1
2
coembed$blacklist_region_fragments[is.na(coembed$blacklist_region_fragments)] <- 0
FeaturePlot(coembed, features = "blacklist_region_fragments", max.cutoff = 500)

黑名单区


Seurat这种基于基因活跃度得分进行细胞类型预测,是否靠谱,开发者提供了如下几个证据

  • 总体预测得分(pbmc.atac$prediction.score.max)高,意味着用scRNA-seq定义细胞类型比较可靠
  • 我们可以在scATC-seq降维结果中
  • 利用相同锚点的贡嵌入分析,发现两类形态能很好的混合
  • 将ATAC-seq数据根据聚类结果构建pseduo bulk, 发现和真实的bulk数据近似

参考资料:

为ChIPseeker增加了一个subset函数

有一天我的师弟提了一个需求:

对于ChIP的下游分析,我很喜欢DiffBind做差异分析然后用ChIPseeker做注释这一套流程(因为ChIPseeker的输入格式是GRange格式,而DiffBind的dba.report输出也恰好是GRange格式,两者可以无缝衔接)。我之前的套路,一直都是无脑输出所有的DiffBind结果,然后放入ChIPSeeker里面做注释,完了之后转成data.frame,然后对data.frame做subset,做后续的GO分析()。但今天突然想对ChIPseeker的注释结果直接做subset,然后分别对上下调的结果画plotAnnoPie等ChipSeeker内置的画图函数。但发现subset无法对ChipSeeker annotatePeak函数的输出结果做筛选。

当然,对于DiffBind的结果用subset,分别提取上下调,然后再注释也是可以的,不过这样感觉更麻烦。

面对需求,如果无法解决提出需求的人,那么就只能写代码实现这个需求了。于是经过一个下午的努力,我给ChIPseeker加上了subset的功能。

先来介绍下如何使用。我们需要用devtools安装最新的ChIPseeker

1
2
3
4
5
BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")
install.packages("ggupset")
install.packages("ggimage")
devtools::install_github("YuLab-SMU/ChIPseeker")
library(ChIPseeker)

我们以测试数据集来演示

1
2
3
4
5
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
peakfile <- system.file("extdata", "sample_peaks.txt", package="ChIPseeker")
peakAnno <- annotatePeak(peakfile, tssRegion=c(-3000, 3000), TxDb=txdb)
peakAnno

此时会输出peakAnno里的描述统计信息

1
2
3
4
5
6
7
8
9
10
11
12
13
14
Annotated peaks generated by ChIPseeker
150/150 peaks were annotated
Genomic Annotation Summary:
Feature Frequency
8 Promoter (<=1kb) 4.666667
9 Promoter (1-2kb) 4.000000
10 Promoter (2-3kb) 2.000000
3 5' UTR 1.333333
2 3' UTR 2.666667
6 Other Exon 8.000000
1 1st Intron 13.333333
7 Other Intron 32.000000
5 Downstream (<=300) 1.333333
4 Distal Intergenic 30.666667

peakAnno可以直接用来画图

1
2
plotAnnoPie(peakAnno)
upsetplot(peakAnno)

Pie

upset

假如你突发奇想,我能不能就看看某一条染色体的作图结果,或者对于MACS2的结果,想根据FDR筛选下peak,然后看下peak的变化。

在之前我们无法直接对peakAnno背后的csAnno对象进行操作,所以需要先对输入的GRanges进行过滤然后再用annotatePeak进行注释,然后画图。

但是现在我们有了subset, 这一切就变得稍微简单起来了

我们可以按照染色体编号进行筛选

1
2
peakAnno_subset <- subset(peakAnno, seqnames == "chr1")
upsetplot(peakAnno_subset, vennpie = TRUE)

可以按照peak宽度进行筛选

1
2
peakAnno_subset <- subset(peakAnno, length > 500)
upsetplot(peakAnno_subset, vennpie = TRUE)

filter

如果筛选之后没有任何结果,那么就会报错

1
2
3
> subset(peakAnno, FDR... < .05)
Error in names(x) <- value :
'names' attribute [2] must be the same length as the vector [1]

那么问题就是,这些筛选条件的名字怎么获知。

只要将peakAnno转成GRanges格式,所以可以看到的列都是你可以筛选的标准。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
gr <- as.GRanges(peakAnno)
head(gr)
GRanges object with 2 ranges and 15 metadata columns:
seqnames ranges strand | length summit tags
<Rle> <IRanges> <Rle> | <integer> <integer> <integer>
[1] chr10 105137981-105138593 * | 614 174 7
[2] chr10 42644417-42645383 * | 968 669 27
X.10.log10.pvalue. fold_enrichment FDR...
<numeric> <numeric> <numeric>
[1] 52.8 15.32 <NA>
[2] 77.15 9.13 0.79
annotation geneChr geneStart geneEnd
<character> <integer> <integer> <integer>
[1] Exon (uc001kwv.3/6877, exon 3 of 11) 10 105127724 105148822
[2] Distal Intergenic 10 42827314 42863493
geneLength geneStrand geneId transcriptId distanceToTSS
<integer> <integer> <character> <character> <numeric>
[1] 21099 1 6877 uc010qqq.2 10257
[2] 36180 2 441666 uc010qey.2 218110

给ChIPseeker增加subset的功能,应该是我第一次在GitHub进行pull requests操作。写这个函数,需要先去理解csAnno对象到底是什么,以及如何保证csAnno这个对象在操作前后不会发生变化。

其实这个函数挺简单的,就是下面几行

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
##' @importFrom S4Vectors subset
##' @importFrom BiocGenerics start
##' @importFrom BiocGenerics end
##' @method subset csAnno
##' @export
subset.csAnno <- function(x, ... ){

index <- paste(seqnames(x@anno),start(x@anno),end(x@anno), sep = "_")
# subset the GRanges
x@anno <- subset(x@anno, ...)
index2 <- paste(seqnames(x@anno),start(x@anno),end(x@anno), sep = "_")

# the tssRgion, level, hsaGenomicAnnotation keep unchanged

# change the detailGenomicAnnotation
x@detailGenomicAnnotation <- x@detailGenomicAnnotation[index %in% index2,]

# change the annotation stat
x@annoStat <- getGenomicAnnoStat(x@anno)

# change peak number
x@peakNum <- length(x@anno)

return(x)

}

因为有很多功能是ChIPseeker已经有了的,比如说getGenomicAnnoStat, 还有一些是S4Vectors和BiocGenerics包提供的,比如说subset, start和end。

三代组装软件miniasm笔记

我们用来练手的文章发表在 Nature Communication ,”High contiguity Arabidopsis thaliana genome assembly with a single nanopore flow cell”, 非常不要脸的说,这篇文章是我师爷实验室发的。

简单讲讲故事内容,就是他们实验室买了一台nanopore仪器,就是下面这台, 目前仪器价格国内是8K左右,当然测序的价格就另说了。如同买台PS4主机,还要买游戏,买个单反,你还得买镜头。仪器只是败家的开始!

nanopore

他们认为三代测序目前有两大问题,测的还不够长以及不够准。nanopore解决了其中一个问题,不够长。Arabidopsis thaliana 当年用一代测序,虽然可以认为是组装的金标准了,但是还是有很多区域是BAC连BAC文库搞不定的,所以就用这台仪器把 Arabidopsis thaliana 测了一波。显然就测一个nanopore,还是已知序列的物种是不可能发文章的,于是他们又用Pacbio sequel测了一波。最后用bionano 光学图谱验证了一次(请大家自行计算要多少钱)。

光测序不行,还得组装对吧。传统的组装方法是想办法利用高深度和随机错误进行纠错,然后用纠错后的长序列进行组装,最后用二代进行纠错。对于一台不错的服务器(20W起步吧)大约花个十天半个月就行。作者或许认为买一台20多w的外设配合不到1w的测序仪可能是太蠢了,于是他用了比较Li Heng大神开发的工具,Minimap+miniasm进行组装,然后用racon+pillon进行纠错,用了一台Macbook Pro 15.6寸花了4天就搞定了,并且和常规工具比较,还算过得去哦。

下面就是正式的分析:

根据文章提供的项目编号”PRJEB21270”, 在European Nucleotide Archive上找到下载地址。

ENA搜索

进入这个页面之后,就可以去下载作者用到的所有数据,我们下载Sequel和MinIon和Illuminia的数据就好了,数据量加起来差不多30G。

1
2
3
4
5
## Sequal
wget -c -q ftp://ftp.sra.ebi.ac.uk/vol1/ERA111/ERA1116568/bam/pb.bam
wget -c -q ftp://ftp.sra.ebi.ac.uk/vol1/ERA111/ERA1116568/bam/pb.bam.bai
## MinION
wget -c -q ftp://ftp.sra.ebi.ac.uk/vol1/ERA111/ERA1116595/fastq/ont.fq.gz

对于Illumina的二代测序,需要用prefetch进行下载

1
2
3
# Illuminia MiSeq
prefetch ERR2173372
fasterq-dump -O . ERR2173372

拿到数据之后,我们就可以用作者提供的分析流程进行重复了。地址为https://github.com/fbemm/onefc-oneasm/wiki/Assembly-Generation

这就是大神的自信,把代码都给你,反正你也看不懂。当然我在重复的时候用的都是最新的软件,所以会有所不同

第一步:拿着80%~90%正确率的原始数据相互比对, 找序列之间的Overlap。这一步,我花了30分钟

1
time ~/opt/biosoft/minimap2/minimap2 -t 10 -x ava-ont ont.fq ont.fq > gzip -1 ont.paf.gz &

第二步:找到Overlap,就能够进行组装了。这一步我花了2分钟

1
2
time ~/opt/biosoft/miniasm/miniasm -f ont.fq ont.paf > ONTmin.gfa &
awk '/^S/{print ">"$2"\n"$3}' ONTmin.gfa | seqkit seq > ONTmin_IT0.fasta &

第三步: 原始的组装结果充满了错误,所以需要进行纠错。纠错分为两种,一种是用三代自身数据,一种是用二代数据进行纠错。当然这两步都是需要的

首先使用三代数据进行纠错,古语有云“事不过三”一般迭代个三次就差不多。这三步,差不多用了1个小时。

1
2
3
4
5
6
7
8
9
# Iteration 1
~/opt/biosoft/minimap2/minimap2 ONTmin_IT0.fasta ont.fq > ONTmin_IT0.paf &
time ~/opt/biosoft/racon/build/bin/racon -t 10 ont.fq ONTmin_IT0.paf ONTmin_IT0.fasta > ONTmin_IT1.fasta &
# Iteration 2
~/opt/biosoft/minimap2/minimap2 ONTmin_IT1.fasta ont.fq > ONTmin_IT1.paf
time ~/opt/biosoft/racon/build/bin/racon -t 10 ont.fq ONTmin_IT1.paf ONTmin_IT1.fasta> ONTmin_IT2.fasta
# Iteration 3
~/opt/biosoft/minimap2/minimap2 ONTmin_IT2.fasta ont.fq > ONTmin_IT2.paf
time ~/opt/biosoft/racon/build/bin/racon -t 10 ont.fq ONTmin_IT2.paf ONTmin_IT2.fasta > ONTmin_IT3.fasta

之后使用二代数据进行纠错。二代数据虽然短,但是测序质量高,所以一般都要用它进行纠错。推荐用30X PCR free的illuminia 测序数据。

Step 1: 数据预处理,过滤低质量短读,去接头。工具很多,常用的是trimmomatic,cutadapter. 我安利一个国内海普洛斯搞的一个工具fastp。

1
2
# data clean
fastp -q 30 -5 -l 100 -i ERR2173372_1.fastq -I ERR2173372_2.fastq -o i1_clean_1.fq -O i1_clean_2.fq

这里标准为:平均质量高于Q30,对5‘端进行低质量碱基删除,保留大于100bp的短读

Step2: 比对,这一步基本都只用了bwa了

1
2
3
# align
bwa index ONTmin_IT3.fasta
bwa mem -t 8 ONTmin_IT3.fasta il_clean_1.fastq il_clean_2.fastq | samtools sort -@ 8 > ONTmin_IT3.bam

step3: 使用比对后的BAM文件进行纠错

1
2
# short read consensus call
java -Xmx16G -jar pilon-1.22.jar --genome ONTmin_IT3.fasta --frags ONTmin_IT3.bam --fix snps --output ONTmin_IT4

二代纠错的时间明显比之前的久,需要一天时间。

大家拿出自己的笔记本实际感受下呗

参考文献

  • nanopore组装拟南芥: High contiguity Arabidopsis thaliana genome assembly with a single nanopore flow cell
  • 不纠错组装: Minimap and miniasm: fast mapping and de novo assembly for noisy long sequences
  • 三代组装软件评测: Comprehensive evaluation of non-hybrid genome assembly tools for third-generation PacBio long-read sequence data

使用purge_haplogs处理基因组杂合区域

FALCON和Canu的组装后会得到一个单倍型融合的基因组,用来表示二倍体基因组。之后,FALCON Unzip和Supernova这类软件进一步处理其中等位基因区域,将这部分区间进行拆分。

当基因组某些区域可能有着比较高的杂合度,这会导致基因组该区域的两个单倍型被分别组装成primary contig, 而不是一个为primary contig, 另一个是associated haplotig. 如果下游分析主要关注于单倍型,这就会导致一些问题。

那么有没有解决方案呢?其实也很好办,就是找到相似度很高的contig,将他们拆分。目前已经有一些软件可以完成类似任务,如 HaploMerger2, Redundans, 这不过这些软件主要处理二代组装结果。

purge_haplogs则是最近开发,用于三代组装的基因组。它根据minimap2的比对结果,通过分析比对read的覆盖度决定谁去谁留。该工具适用于单倍型组装软件,例如 Canu, FALCON或 FALCON-Unzip primary contigs, 或者是分相后的二倍体组装(Falcon-Unzip primary contigs + haplotigs 。

它的工作流程如下图所示。一般只需要两个输入文件,组装草图(FASTA格式) 和 比对的BAM文件。同时还可以提供重复序列注释的BED文件,辅助处理高重复的contig。

分析流程

建议: 用原来用于组装的read进行比对。对于多个匹配的read,建议采取random best,也就是随便找一个。

软件安装

purge_haplotigs依赖软件比较多,手动安装会很麻烦,但是他可以直接用bioconda装

1
2
3
conda create -n purge_haplotigs_env
conda activate purge_haplotigs_env
conda install purge_haplotigs

安装完成后需要一步测试

1
purge_haplotigs test

简明教程

数据准备。 需要下载的数据集分为两个部分,一个是FALCON-Unzip后的primary contig 和 halplotigs. 另一个则是已经比完后的BAM文件

1
2
3
4
5
6
7
mkdir purge_haplotigs_tutorial
cd purge_haplotigs_tutorial
wget https://zenodo.org/record/841398/files/cns_h_ctg.fasta
wget https://zenodo.org/record/841398/files/cns_p_ctg.aligned.sd.bam # 1.7G
wget https://zenodo.org/record/841398/files/cns_p_ctg.aligned.sd.bam.bai
wget https://zenodo.org/record/841398/files/cns_p_ctg.fasta
wget https://zenodo.org/record/841398/files/cns_p_ctg.fasta.fai

当然我们不可能直接就拿到比对好的BAM文件,我们一般是有组装后的基因组以及用于组装的subread,假设这两个文件命名为, genome.fa 和 subreads.fasta.gz.

官方提供的新比对代码

1
2
minimap2 -t 4 -ax map-pb genome.fa subreads.fasta.gz --secondary=no \
| samtools sort -m 1G -o aligned.bam -T tmp.ali

如下是旧版代码

1
2
3
minimap2 -ax map-pb genome.fa subreads.fasta.gz \
| samtools view -hF 256 - \
| samtools sort -@ 8 -m 1G -o aligned.bam -T tmp.ali

如果你有二代测序数据,也可以用BWA-MEM进行比对得到BAM文件。

第一步:使用purge_haplotigs readhist从BAM中统计read深度,绘制柱状图。

1
2
3
4
5
# 新
purge_haplotigs hist -b aligned.bam -g genome.fasta -t 20
# 旧
# purge_haplotigs readhist -b aligned.bam -g genome.fasta -t 20
# -t 线程数, 不宜过高,过高反而没有效果。

也就是下图,你能明显的看到图中有两个峰,一个是单倍型的覆盖度,另一个二倍型的覆盖度,

高杂合基因组read-depth histogram

你可能还想知道高纯合基因组是什么样的效果,我也找了一个纯合的物种做了也做了read-depth 柱状图,

纯合基因组read-depth histogram

之后你需要根据read-depth 柱状图 确定这两个峰的位置用于下一步。下面是两个例子。对于我们则是,20,65,190.

两个例子

第二步: 根据read-depth信息选择阈值。

1
2
3
4
# 新
purge_haplotigs contigcov -i cns_p_ctg.aligned.sd.bam.gencov -o coverage_stats.csv -l 20 -m 75 -h 190
# 旧
# purge_haplotigs contigcov -i cns_p_ctg.aligned.sd.bam.gencov -o coverage_stats.csv -l 20 -m 75 -h 190

这一步生成的文件是”coverage_stats.csv”

第三步:区分haplotigs.

1
purge_haplotigs purge  -g cns_p_ctg.fasta  -c coverage_stats.csv  -b cns_p_ctg.aligned.sd.bam  -t 4  -a 60

这一步会得到如下文件

  • curated.artefacts.fasta:无用的contig,也就是没有足够覆盖度的contig.
  • curated.fasta:新的单倍型组装
  • curated.haplotigs.fasta:从原本组装分出来的haplotigs
  • curated.reassignments.tsv: 单倍型的分配信息
  • curated.contig_associations.log: 运行日志, 下面是其中一个记录,表示000004F_004和000004F_027是000004F_017的HAPLOTIG, 而000004F_017和000004F_013又是000004F,的HAPLOTIG。
1
2
3
4
000004F,PRIMARY -> 000004F_013,HAPLOTIG
-> 000004F_017,HAPLOTIG
-> 000004F_004,HAPLOTIG
-> 000004F_027,HAPLOTIG

在”curated.reassignments.tsv”文件中有6列

  • reassigned_contig: 用于比较的contig
  • top_hit_contig: 最好的被比对的contig
  • second_hit_contig: 第二个被比对的contig
  • best_match_coverage: 最好的匹配覆盖度
  • max_match_coverage : 最高的匹配深度
  • reassignment: 标记为haplotype 还是 repeat,或者是keep

由于我们用的是单倍型组装primary contigs而不是二倍体组装的parimary + haplotigs, 因此我们需要将FALCON_Unzip的haplotgi合并到重新分配的haplotigs中,这样子我们依旧拥有二倍体组装结果

1
cat cns_h_ctg.fasta >> curated.haplotigs.fasta

检查dotplots

如果在第三步purge_haplotigs purge中添加了-d/--dotplots参数,即为每个reassigned_contigs和unassigned_contigs生成用于人工检查的共线性图,那么在最终的输出结果中会有两个文件目录

  • dotplots_reassigned_contigs
  • dotplots_unassigned_contigs

官方提供了5个可能出现的共线性图

第一种: Haplotig,最佳情况,完美的共线性关系

Haplotig

第二种: 大部分是haplotigs, 说明这个contig部分是二倍型,部分是单倍型,可能是半合子(hemizygous)

mostly haplotig

第三种: Haplotigs里有大量的串联重复

tandem repeat

第四种: Haplotigs是回文序列(palindrome)

palindrome

第五种: contig从string graph中knots产生,这种情况不算是haplotigs,但是对于短读序列的比对会造成麻烦

knots

你可能需要看大量的图才能有感觉,到底应该把哪些Purge_haplotigs错误认为是haplotig的contig放回到primary contig中。

原理介绍

为什么第一步的 read 覆盖深度分析能判断基因组是否冗余呢?这是因为对于坍缩的单倍型,那么含有等位的基因的read只能比对到该位置上,而如果杂合度太高被拆分成两个不同的contig,那么含有等位的基因的read就会分别比对到不同的read上,导致深度降低一半。下图A就是一个典型的包含冗余基因组的read覆盖度分布

read-深度分析

分析流程的第二步的任务就是人工划分出如下图B部分,绿色的部分是坍缩单倍型contig,蓝色的部分是潜在的冗余contig。之后,分析流程会计算这些区域中的contig的覆盖度。 对于绿色部分中的contig,如果覆盖度低于80%, 会进行标记用于后续分析。如果深度非常低,那么很有可能就是组装引入错误,深度非常高的部分基本就是重复序列或者是细胞器的contig,这些黄色的contig可以在后续的组装出分开。

划分区间

第三步就是同源序列进行识别和分配。所有标记的contig之后会用Minimap2在整个组装进行搜索,寻找相似度较高的离散区间(如下图C)。如果一个Contig的联配得分大于阈值(默认70%), 那么就会被标记为haplotigs. 如果一个contig的最大联配得分大于阈值(默认250%), 会被标记成重复序列,这有可能是潜在的有问题contig,或许是坍缩的contig或者低复杂度序列。

移除haplotigs

推荐阅读

使用TransDecoder寻找转录本中的编码区

TransDecoder能够从转录本序列中鉴定候选编码区。这些转录本序列可以来自于Trinity的从头组装,或者来自于Cufflinks或者StringTie的组装结果。

软件安装

https://github.com/TransDecoder/TransDecoder/releases下载最新版的TransDecoder,以v5.5.0为例

1
2
3
4
mkdir -p ~/opt/biosoft && cd ~/opt/biosoft
wget https://github.com/TransDecoder/TransDecoder/archive/TransDecoder-v5.5.0.zip
unzip TransDecoder-v5.5.0.zip
mv TransDecoder-TransDecoder-v5.5.0 TransDecoder-v5.5.0

运行TransDecoder

我们从cufflinks或stringtie输出的gtf文件开始分析流程,因此你会有两个输入文件

  • transcripts.gtf: 记录预测转录本的GTF文件
  • genome.fasta: 参考基因组序列

第一步: 从GTF文件中提取FASTA序列

1
~/opt/biosoft/TransDecoder-v5.5.0/util/gtf_genome_to_cdna_fasta.pl transcripts.gtf genome.fasta > transcripts.fasta

第二步: 将GTF文件转成GFF3格式

1
~/opt/biosoft/TransDecoder-v5.5.0/util/gtf_to_alignment_gff3.pl transcripts.gtf > transcripts.gff3

第三步: 预测转录本中长的开放阅读框, 默认是100个氨基酸,可以用-m修改

1
~/opt/biosoft/TransDecoder-v5.5.0/TransDecoder.LongOrfs -t transcripts.fasta

第四步: 使用DIAMOND对上一步输出的transcripts.fasta.transdecoder.pep在蛋白数据库中进行搜索,寻找同源证据支持

1
2
3
4
5
6
7
# 下载数据并解压缩
wget ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz
gunzip uniprot_sprot.fasta.gz
# 建立索引
diamond makedb --in uniprot_sprot.fasta --db uniprot_sprot.fasta
# BLASTP比对
diamond blastp -d uniprot_sprot.fasta -q transcripts.fasta.transdecoder_dir/longest_orfs.pep --evalue 1e-5 --max-target-seqs 1 > blastp.outfmt6

关于DIAMOND的使用,参考这篇说明DIAMOND: 超快的蛋白序列比对软件

第五步: 预测可能的编码区

1
2
3
~/opt/biosoft/TransDecoder-v5.5.0/TransDecoder.Predict \
-t transcripts.fasta \
--retain_blastp_hits blastp.outfmt6

第六步: 生成基于参考基因组的编码区注释文件

1
2
3
4
~/opt/biosoft/TransDecoder-v5.5.0/util/cdna_alignment_orf_to_genome_orf.pl \
transcripts.fasta.transdecoder.gff3 \
transcripts.gff3 \
transcripts.fasta > transcripts.fasta.transdecoder.genome.gff3

最终输出文件如下:

  • transcripts.fasta.transdecoder.pep: 最终预测的CDS对应的蛋白序列
  • transcripts.fasta.transdecoder.cds: 最终预测的CDS序列
  • transcripts.fasta.transdecoder.gff3: 最终ORF对应的GFF3
  • transcripts.fasta.transdecoder.bed: 以BED格式存放ORF位置信息
  • transcripts.fasta.transdecoder.genome.gff3: 基于参考基因组的GFF3文件

其中BED和GFF3可以放到IGV上展示,手动检查下结果

假如是Trinity从头预测的转录本,没有参考基因组,那么就运行第三步,第四步和第五步

在transcripts.fasta.transdecoder.cds文件中,每个fasta序列的头信息部分会有一个ORF type,分为如下几个类型

  • ‘complete’: 包含起始密码子和终止密码子
  • ‘5prime_partial’: 缺失起始密码子, 可能只有部分的N端序列
  • ‘3prime_partial’: 缺失终止密码子, 可能只有部分的C端序列
  • ‘internal’: 意味着同时是5prime-partial和3prime-partial

通常而言,我们会用complete的cds寻找同源证据,然后选择高可信度的序列用于训练,而不是哪些特别长且没有已知蛋白支持的序列。

参考资料

DIAMOND:超快的蛋白序列比对软件

今天用BLASTX将我的转录本序列在UniProt蛋白数据库(700w条序列)中搜索,80个线程,过了1小时大概就分析1000条吧。实在是有点慢,于是我想到之前耳闻的DIAMOND,据说速度非常快,于是我测试了下。没想到,这工具居然那么快。

根据DIAMOND介绍,它有以下特点

  • 比BLAST快500到20,000倍
  • 长序列的移框联配分析(frameshift alignment)
  • 资源消耗小,普通台式机和笔记本都能运行
  • 输出格式多样

我就看中它一点,速度快。

软件安装异常的简单,因为提供了预编译的64位可执行文件

1
2
3
4
5
6
wget http://github.com/bbuchfink/diamond/releases/download/v0.9.25/diamond-linux64.tar.gz
tar xzf diamond-linux64.tar.gz
# 有root全新啊
sudo mv diamond /usr/local/bin
# 无root权限, ~/bin是自己当前目录下
mv diamond ~/bin

因为 diamon的功能就是将蛋白或者翻译后的核苷酸和蛋白数据库进行比对,没有BLAST那么多功能,所以软件使用也是异常的简单。

第一步: 先从NCBI上下载蛋白数据库。 NR库是NCBI的非冗余蛋白数据库,

1
2
wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz
gunzip nr.gz

也可以从ftp://ftp.ncbi.nlm.nih.gov/refseq/release/plant/下载植物的蛋白数据库

第二步: 建库。就两个参数,--in nr输入文件,--db nr 输出的数据库前缀. 氨基酸序列中的结尾可以有”*”

1
diamond makedb --in nr --db nr

: 假如要根据GFF提取蛋白序列,一定要注意输出的氨基酸序列中不能有”.”在序列中,否则会报错。可以通过seqkit grep -s -vrp '"\."' input.fa > output.fa 进行过滤。

第三步: 搜索。就两个子命令,blastp和blastx,前者比对蛋白,后者比对DNA序列

1
2
diamond blastx --db nr -q reads.fna -o dna_matches_fmt6.txt
diamond blastp --db nr -q reads.faa -o protein_matches_fmt6.txt

-q/--query输入检索序列,--out/-o输出文件,默认以--outfmt 6输出结果和BLAST+的--outfmt 6结果一致。

注意事项:

  • 默认参数主要是针对短序列,对于比较长的序列,使用--sensitive--more-senstive提高敏感度。
  • 默认的e-value阈值是0.001, 而BLAST是10,因此会比BLAST结果更加严格

性能优化:

  • 设置比较低的-e参数
  • 设置-k参数,减少输出的联配数目。这会降低临时文件大小和最终结果
  • --top会输出得分比最好的分数低一定百分比的结果,
  • --compress 1: 输出结果会以gzip进行压缩

参考文献

Benjamin Buchfink, Chao Xie, and Daniel H. Huson. Fast and sensitive protein alignment
using diamond. Nature methods, 12(1):59–60, Jan 2015.

使用MAKER进行全基因组基因注释-基础篇

maker

在基因组注释上,MAKER算是一个很强大的分析流程。能够识别重复序列,将EST和蛋白序列比对到基因组,进行从头预测,并在最后整合这三个结果保证结果的可靠性。此外,MAKER还可以不断训练,最初的输出结果可以继续用作输入训练基因预测的算法,从而获取更高质量的基因模型。

Maker的使用比较简单,在软件安装成后,会有一个”data”文件夹存放测试数据

1
2
ls ~/opt/biosoft/maker/data
dpp_contig.fasta dpp_est.fasta dpp_protein.fasta hsap_contig.fasta hsap_est.fasta hsap_protein.fasta te_proteins.fasta

以”dpp”开头的数据集为例,protein表示是同源物种的蛋白序列,est是表达序列标签,存放的是片段化的cDNA序列,而contig则是需要被预测的基因组序列。

让我们新建一个文件夹,并将这些测试数据拷贝过来。

1
2
mkdir test01 ; cd test01
cp ~/opt/biosoft/maker/data/dpp* .

由于基因组注释设计到多个程序,多个步骤,每个步骤可能都有很多参数需要调整,因此就需要建立专门的配置文件用来告诉maker应该如何控制流程的运行。

如下步骤创建三个以ctl结尾的配置文件

1
2
3
~/opt/biosoft/maker/bin/maker -CTL
ls *.ctl
maker_bopts.ctl maker_exe.ctl maker_opts.ctl
  • maker_exe.ctl: 执行程序的路径
  • maker_bopt.ctl: BLAST和Exonerat的过滤参数
  • maker_opt.ctl: 其他信息,例如输入基因组文件

maker_exe.ctl和maker_bopt.ctl可以简单用less查看,可不做修改,maker_opt.ctl是主要调整的对象。 使用vim maker_opt.ctl修改如下内容

1
2
3
4
genome=dpp_contig.fasta
est=dpp_est.fasta
protein=dpp_protein.fasta
est2genome=1

修改完之后多花几分钟看看每个参数的设置,尽管很枯燥,但是考虑这个工具你可能会反复多次使用,所以这点时间是一定要花的。

随后就可以在当前路径运行程序

1
~/opt/biosoft/maker/bin/maker &> maker.log &

输出结果见”dpp_contig.maker.output”, 重点是”dpp_contig_master_datastore_index.log”文件,由于maker会拆分数据集并行计算,因此该文件记录总体的运行情况,需要关注其中是否有”FAILED”,”RETRY”,”SKIPPED_SAMLL”,”DIED_SIPPED_PERMANET”,因为这意味着有些数据出于某些原因没有运算。

最后,我们需要将并行运算的结果进行整合,导出GFF文件, 转录本序列和蛋白序列

1
2
~/opt/biosoft/maker/bin/fasta_merge -d dpp_contig_master_datastore_index.log
~/opt/biosoft/maker/bin/gff3_merge -d dpp_contig_master_datastore_index.log

在该目录下就会出现, “dpp_contig.all.gff”, “dpp_contig.all.maker.proteins.fasta”,”dpp_contig.all.maker.transcripts.fasta”

其中GFF文件就需要用IGV,JBrowse, Apollo下展示来检查下注释是否正确。

附录

软件安装:MAKER可以免费用于学术用途,但是未经许可不可商用。目前有两个版本2018年5月4日更新的2.31.10和测试版3.01.02.出于稳定性考虑,安装前者。后续假设已经在http://yandell.topaz.genetics.utah.edu/cgi-bin/maker_license.cgi进行登记,并且下载了压缩包”maker-2.31.10.tgz”

先检查下自己的系统情况,看需要补充哪些库

1
2
3
tar xf maker-2.31.10.tgz
cd maker/src
perl Build.PL

这一步之后会罗列出后续需要运行的命令来完成安装

1
2
3
4
./Build installdeps
./Build installexes
./Build install
./Build status

参考资料

  • Genome Annotation and Curation Using MAKER and MAKER-P

三代组装软件Canu使用

Canu简介

Canu是Celera的继任者,能用于组装PacBio和Nanopore两家公司得到的测序结果。

Canu分为三个步骤,纠错,修整和组装,每一步都差不多是如下几个步骤:

  • 加载read到read数据库,gkpStore
  • 对k-mer进行技术,用于计算序列间的overlap
  • 计算overlap
  • 加载overlap到overlap数据库,OvlStore
  • 根据read和overlap完成特定分析目标
    • read纠错时会从overlap中挑选一致性序列替换原始的噪声read
    • read修整时会使用overlap确定read哪些区域是高质量区域,哪些区域质量较低需要修整。最后保留单个最高质量的序列块
    • 序列组装时根据一致的overlap对序列进行编排(layout), 最后得到contig。

这三步可以分开运行,既可以用Canu纠错后结果作为其他组装软件的输入,也可以将其他软件的纠错结果作为Canu的输入,因此下面分别运行这三步,并介绍重要的参数。

几个全局参数:genomeSize设置预估的基因组大小,这用于让Canu估计测序深度; maxThreads设置运行的最大线程数;rawErrorRate用来设置两个未纠错read之间最大期望差异碱基数;correctedErrorRate则是设置纠错后read之间最大期望差异碱基数,这个参数需要在 组装 时多次调整;minReadLength表示只使用大于阈值的序列,minOverlapLength表示Overlap的最小长度。提高minReadLength可以提高运行速度,增加minOverlapLength可以降低假阳性的overlap。

组装实战

数据下载

数据来自于发表在 Nature Communication 的一篇文章 “High contiguity Arabidopsis thaliana genome assembly with a single nanopore flow cell“。 这篇文章提供了 Arabidopsis thaliana KBS-Mac-74 的30X短片段文库二代测序、PacBio和Nanopore的三代测序以及Bionano测序数据, 由于拟南芥的基因组被认为是植物中的金标准,因此文章提供的数据适合非常适合用于练习。根据文章提供的项目编号”PRJEB21270”, 在European Nucleotide Archive上找到下载地址。

ENA搜索

1
2
3
4
5
6
7
## PacBio Sequal
wget -c -q ftp://ftp.sra.ebi.ac.uk/vol1/ERA111/ERA1116568/bam/pb.bam
## MinION
wget -c -q ftp://ftp.sra.ebi.ac.uk/vol1/ERA111/ERA1116595/fastq/ont.fq.gz
# Illuminia MiSeq
wget -c -q ftp://ftp.sra.ebi.ac.uk/vol1/ERA111/ERA1116569/fastq/il_1.fq.gz
wget -c -q ftp://ftp.sra.ebi.ac.uk/vol1/ERA111/ERA1116569/fastq/il_2.fq.gz

下载的PacBio数据以BAM格式存储,可以通过安装PacBio的smrtlink工具套装,使用其中的bam2fasta工具进行转换.

1
2
3
4
# build index for convert
~/opt/biosoft/smrtlink/smrtcmds/bin/pbindex pb.bam &
# convert bam to fasta
~/opt/biosoft/smrtlink/smrtcmds/bin/bam2fasta -o pb pb.bam &

其实samtools fasta也可以将bam转成fasta文件,并且不影响之后的组装。

PacBio的smrtlink工具套装大小为1.4G,不但下载速度慢,安装也要手动确认各种我不清楚的选项, 唯一好处就是工具很全。

运行Canu

第一步:纠错。三代测序本身错误率高,使得原始数据充满了噪音。这一步就是通过序列之间的相互比较纠错得到高可信度的碱基。主要调整两个参数

  • corOutCoverage: 用于控制多少数据用于纠错。比如说拟南芥是120M基因组,100X测序后得到了12G数据,如果只打算使用最长的6G数据进行纠错,那么参数就要设置为50(120m x 50)。设置一个大于测序深度的数值,例如120,表示使用所有数据。
1
2
3
4
5
6
canu -correct \
-p ath -d pb_ath \
Threads=10 gnuplotTested=true\
genomeSize=120m minReadLength=2000 minOverlapLength=500\
corOutCoverage=120 corMinCoverage=2 \
-pacbio-raw pb.fasta.gz

可以将上述命令保存到shell脚本中进行运行, nohup bash run_canu.sh 2> correct.log &.

注: 1.7版本里会没有默认没有安装gnuplot出错,因此gnuplotTested=true 可以跳过检查。

第二步:修整。这一步的目的是为了获取更高质量的序列,移除可疑区域(如残留的SMRTbell接头).

1
2
3
4
5
canu -trim \
-p ath -d pb_ath
maxThreads=20 gnuplotTested=true\
genomeSize=120m minReadLength=2000 minOverlapLength=500\
-pacbio-corrected ath/pb_ath.correctedReads.fasta.gz

第三步: 组装。在前两步获得高质量的序列后,就可以正式进行组装. 这一步主要调整的就是纠错后的序列的错误率, correctedErrorRate,它会影响utgOvlErrorRate。这一步可以尝试多个参数,因为速度比较块。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# error rate 0.035
canu -assemble \
-p ath -d ath-erate-0.035 \
maxThreads=20 gnuplotTested=true \
genomeSize=120m\
correctedErrorRate=0.035 \
-pacbio-corrected atg/pb_ath.trimmedReads.fasta.gz
# error rate 0.050
canu -assemble \
-p ath -d ath-erate-0.050 \
maxThreads=20 gnuplotTested=true \
genomeSize=120m\
correctedErrorRate=0.050 \
-pacbio-corrected atg/pb_ath.trimmedReads.fasta.gz

最后输出文件下的ath.contigs.fasta就是结果文件。

一些宝贵的建议

Nanopore组装

对于Nanopore数据,使用Canu组装并不是一个非常好的选择,我曾经以一个600多M的物种100X数据进行组装,在120线程,花了整整一个多月的时间,尽管它的组装效果真的是很好。

  • rawErrorRate: 从0.144 调整到 0.12 或者更低,速度会提高5到10倍
  • readSamplingCoverage, readSamplingBias: 可以抽取部分数据进行纠错,而不是全部数据
  • corOutCoverage=30: 默认的30X或者40X的组装效果其实不错
  • -fast: 对于1Gb以下的物种,可以加上该参数,会明显提高速度

高杂合物种组装

对于高杂合物种的组装,Canu建议是用 batOptions=-dg 3 -db 3 -dr 1 -ca 500 -cp 50参数尽量分出两套单倍型,然后对基因组去冗余。

batOptions表示传递后续的参数给组装软件bogart, -dg 3 -db3降低自动确定阈值时的错误率离差(deviation),从而更好的分开单倍型。-dr 1 -ca 500 -cp 50会影响错误组装的拆分,对于一个模棱两可的contig,如果至少另一条可选路径的overlap长度至少时500bp,或者说另一条可选路径时在长度上和当前最佳路径存在50%的差异,那么就将contig进行拆分。

关于杂合物种组装的讨论,参考https://github.com/marbl/canu/issues/201#issuecomment-233750764

购买SSD避免服务器IO瓶颈

如果你的服务器线程数很多,你在普通的机械硬盘上运行组装,而且你的系统还是CentOS,那么你需要调整一个参数,避免其中一步的IO严重影响服务器性能。

Canu通过两个策略进行并行,bucketizing (‘ovb’ 任务) 和 sorting (‘ovs’ 任务)。 bucketizing会从1-overlap读取输出的overlap,将他们复制一份作为中间文件。sorting一步将这些文件加载到内存中进行排序然后写出到硬盘上。 如果你的overlap输出特别多,那么该步骤将会瞬间挤爆的你的IO.

为了避免悲剧发生,请增加如下参数: ovsMemory=16G ovbConcurrency=15 ovsConcurrency=15, 也就是降低这两步同时投递的任务数,缓解IO压力。

如何用不同服务器处理同一个任务

overlap这一步时间久,如果并没有服务器集群,而是有多台服务器,可以参考如下方法进行数据并行处理(必须要安装相同的Canu版本)

假如Canu任务中的prefix参数为coc, 用scp进行数据的传递

1
2
3
4
scp -r coc.seqStore wangjw@10.10.87.132:/data1/wjw/Coc/Canu
scp -r unitigging/1-overlapper/overlap.sh wangjw@10.10.87.132:/data1/wjw/Coc/Canu
scp -r unitigging/0-mercounts/ wangjw@10.10.87.132:/data1/wjw/Coc/Canu

到目标服务器的canu目录下

1
2
3
4
mkdir unitigging
mkdir -p unitigging/1-overlapper
mv 0-mercounts/ unitigging/
mv overlap.sh unitigging/1-overlapper/

修改overlap.sh里的bin路径,指定当前服务器的canu路径

之后运行任务即可

1
2
cd unitigging/1-overlapper/
./overlap.sh 77 &> 77.log

运行完任务之后,将目前服务器unitigging/1-overlapper/001目录下的000077.oc 000077.ovb 000077.stats文件传送回原来的服务器unitigging/1-overlapper/001下即可。

参考资料