相关文章推荐
玩足球的闹钟  ·  脑肿瘤有哪些常见症状? | ...·  1 月前    · 
不敢表白的啤酒  ·  郭敬明回应"日韩"微博 ...·  5 月前    · 
大气的针织衫  ·  物理与天文学院原子分子光物理（AMO）与量子 ...·  10 月前    · 
高大的金针菇  ·  USDT稳赚25亿美金PayPal盯上的稳定 ...·  1 年前    · 
飞翔的枕头  ·  看剧吧-《乐高大电影2原声版》喜剧片高清完整 ...·  1 年前    · 
小百科  ›  洲更的第二大脑
伤情的香菜
9 月前

第一步: 在 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数据组装的软件，目前该工具尚未发表，除了 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.fastq 或 xxx.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
  • nanopolish
  • racon

    • 由于拟南芥的基因组比较小，我分别用了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的结果

    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
  • https://github.com/xiaochuanle/NECAT
  • 从ATAC-seq中估计RNA-seq表达水平，即从ATAC-seq reads定量基因表达活跃度
  • 使用LSI学习ATAC-seq数据的内部结构
  • 鉴定ATAC-seq和RNA-seq数据集的锚点
  • 数据集间进行转移，包括聚类的标签，在ATAC-seq数据中推测RNA水平用于共嵌入分析

    • 数据下载

    测试数据下载地址:

  • scATAC-seq : h5格式

  • scATAC-seq metadata : csv文件

  • scRNA-seq : h5格式

    可以复制下载链接到浏览器下载，也可以直接在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))

    从上面的结果中，你可能会发现某些细胞只有在一类技术中存在。举个例子，从巨噬细胞(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数据近似

    • 参考资料：

     
    推荐文章
    玩足球的闹钟  ·  脑肿瘤有哪些常见症状? | Dana-Farber Cancer Institute
    1 月前
    不敢表白的啤酒  ·  郭敬明回应&quot;日韩&quot;微博 郭采洁调侃&quot;小四&quot;穿高跟鞋--文化--人民网
    5 月前
    大气的针织衫  ·  物理与天文学院原子分子光物理（AMO）与量子精密测量方向博士后 ...
    10 月前
    高大的金针菇  ·  USDT稳赚25亿美金PayPal盯上的稳定币是门好生意_手机新浪网
    1 年前
    飞翔的枕头  ·  看剧吧-《乐高大电影2原声版》喜剧片高清完整版免费在线观看
    1 年前
    Link管理   ·   51好读   ·   Sov5搜索   ·   小百科
    小百科 - 百科知识指南