2019-08-23

Linux

如何安装perl模块

由于生物信息早期最多用的语言是perl，因此不可避免就要用别人的perl脚本或者基于perl的项目来处理数据。

使用perl脚本和使用其他编程语言的脚本没啥不同，毕竟你只要传入参数，它就能给你结果。因此对于我们这些不用perl写脚本，只需要调用的人而言，唯一要学会的事情就是**如何安装perl的模块”。

关于perl模块安装，最古老的方法就是使用perl -MCPAN -e shell或者是cpan（两者等价），这也是我最先接触的方法，这里介绍如何使用local::lib和cpanm实现非root权限安装perl模块。

使用系统自带的perl

安装任何软件最怕遇到的问题就是权限问题，因此我们需要先安装local::lib，使得我们能够将perl模块安装到任何地方，简单的说就是安装到我们的家目录下

第一步，下载源代码进行编译安装

wget https://cpan.metacpan.org/authors/id/H/HA/HAARG/local-lib-2.000024.tar.gz
tar xf local-lib-2.000024.tar.gz
cd local-lib-2.000024
perl Makefile.PL --bootstrap=~/opt
make test && make install

第二步：设置环境变量，使得perl在安装模块的时候会优先使用我们指定的路径

1	echo 'eval "$(perl -I$HOME/opt/lib/perl5 -Mlocal::lib=$HOME/opt)"' >> ~/.bashrc

先用perl -I$HOME/opt/lib/perl5 -Mlocal::lib=$HOME/opt表示运行前先添加$HOME/opt/lib/perl5到自己的搜索路径@INC中，然后传入参数$HOME/opt执行模块local::lob，这个模块的执行结果会输出如下内容

Attempting to create directory /home6/wangjw/opt
PATH="/home/zgxu/opt/bin${PATH:+:${PATH}}"; export PATH;
PERL5LIB="/home/zgxu/opt/lib/perl5${PERL5LIB:+:${PERL5LIB}}"; export PERL5LIB;
PERL_LOCAL_LIB_ROOT="/home/zgxu/opt${PERL_LOCAL_LIB_ROOT:+:${PERL_LOCAL_LIB_ROOT}}"; export PERL_LOCAL_LIB_ROOT;
PERL_MB_OPT="--install_base \"/home/zgxu/opt\""; export PERL_MB_OPT;
PERL_MM_OPT="INSTALL_BASE=/home/zgxu/opt"; export PERL_MM_OPT;

这些就作为eval的参数进行执行，也就是说你重启终端后后，PERL5LIB PERL_LOCAL_LIB_ROOT,PERL_MB_OPT,PERL_MM_OPT这几个变量就会重新设置，以此保证你后续安装perl模块时，会优先安装到自己的选择的目录

第三步：安装cpam. 由于之前已经配置了local::lib，因此perl编译的工具都会默认安装到~/opt目录下

wget https://cpan.metacpan.org/authors/id/M/MI/MIYAGAWA/App-cpanminus-1.7043.tar.gz
tar xf App-cpanminus-1.7043.tar.gz
cd App-cpanminus-1.7043
perl Makefile.PL
make test && make install

第四步：使用国内镜像提高下载速度，可以通过别名的方式实现

1	echo 'alias cpanm="cpanm --mirror http://mirrors.163.com/cpan --mirror-only"' >>~/.bashrc

之后便可以使用cpanm Module::Name安装任意的软件了。

自己编译一个perl

自己编译Perl的好处就在于之后的perl模块都会安装到自己的Perl目录下，而不会对系统造成影响。

cd ~/src
wget -4 http://www.cpan.org/src/5.0/perl-5.26.1.tar.gz
tar xf perl-5.26.1.tar.gz
cd perl-5.26.1
./Configure -des -Dprefix=$HOME/opt/sysoft/perl-5.26.1
make test
make install

然后用perl -e '{print "$_\n" foreach @INC}'会发现perl只会在自己的目录~/opt/sysoft/perl-5.26.1下查找模块。那么使用cpanm Module::Name安装的任何包都只会安装到~/opt/sysoft/perl-5.26.1下，你也不需要安装local::lib了

conda的perl和系统的perl冲突

有一次我遇到这个问题

1	perl: symbol lookup error: /home/wangjw/perl5/lib/perl5/x86_64-linux-thread-multi/auto/Cwd/Cwd.so: undefined symbol

这个问题是因为用系统perl安装的软件被conda的perl优先查找到导致，用perl -V和perl -e '{print "$_\n" foreach @INC}'可以发现conda的perl查找路径低于我为系统perl安装的路径，解决方案如下

1	export PERL5LIB=""

2019-08-22

文章重现

3分的水稻转录组分析实战

水稻转录组分析实战

文章标题: Genes, pathways and transcription factors involved in seedling stage chilling stress tolerance in indica rice through RNA-Seq analysis

该文章发表在BMC plant biology，植物科学2区SCI，近四年IF稳定在3分左右。这篇文章的核心部分就是测了好几个转录组，然后做了一点实验验证。

文章简介

本部分内容由我的一位师弟完成

背景

世界50％以上的人口将稻米作为饮食中碳水化合物的主要来源。这是因为水稻的种植范围很广，从山地到平原地区均有种植。正因如此，水稻的一生会面临多种生物胁迫（如病原菌，昆虫等）和非生物胁迫（如寒冷，高温，盐碱等），因此解析水稻的抗逆机制，培育具有多种抗逆性状的水稻具有很大的应用价值。而本文的作者是印度人，印度大部分地区的气候具有明显的旱雨两季，旱季的水稻在苗期易受冷（寒）害的影响。随着RNA-seq手段越来越成熟，在比较不同转录组的差异中发挥了重要作用，这对做基因功能的实验室来讲无疑是巨大的福音。作者就想利用RNA-seq的手段，寻找耐寒品种（Geetanjali）和冷敏感品种（Shahabhagidhan）在转录组水平的差异表达的基因，以期在水稻耐寒过程中鉴定到新颖的转录本，并能获得对基因不同表达水平的认识。

方法

总体思路，对耐寒和冷敏感的品种在冷处理过程中的不同时间点取样，并进行转录组分析，分析完成后，用qRT-PCR来验证表达上调和下调的基因，进一步confirm分析结果。

结果

基于低温胁迫处理的5个不同时间点的比较转录组分析显示：在冷敏感品种中有更多的下调基因，在耐寒品种中有更多上调的基因。在冷敏感品种（CSV）中检测到总计13930个差异基因，而在耐寒品种（CTV）中检测到10599个差异基因，且随时间的延长，差异基因的数目也增加。GO富集分析揭示了18个CSV项和28个CTV项显著参与到molecular function, cellular component and biological process。GO分类揭示了耐寒过程中的差异基因显著参与转录调控，活性氧爆发，脂质结合，催化和水解酶活性的作用。KEGG注释通路揭示了更多数量的基因在耐寒过程中调节不同的通路来发挥作用。

结论

结论真的是简单到吐血，真的是简单到吐血，真的是简单到吐血。就一句话：通过对耐寒品种的转录组分析，我们检测到了涉及到抗逆过程的基因，通路和转录因子。

转录组分析重现

转录组分析环境的搭建和一个分析流程全过程可以哔哩哔哩上搜索”徐洲更”, 即可看我录制的视频。对于学习过程中出现的问题，可以参加我们的付费答疑群，活动见https://mp.weixin.qq.com/s/w47Hc2BmRCLpQMlRZ4rCCw。

本次文章重现涉及到的软件如下:

命令行：prefetch, fasterq-dump, hisat2, samtools, featureCounts, gffread
R：DESeq2, clusterProfiler

1 2	conda create -n rna-seq sra-tools fastqc fastp hisat2 samtools subread gffread multiqc conda activate rna-seq

数据下载

文章用到的RNA-seq数据编号，PRJNA288892，一共有12个样本。我以这个项目编号新建一个文件夹，开始本次的分析

1	mkdir -p PRJNA288892 && cd PRJNA288892

接下来是下载文章用到的数据，编号存在SRR_Acc_List.txt

SRR2089751
SRR2089753
SRR2089756
SRR2093937
SRR2093938
SRR2093939
SRR2093948
SRR2093955
SRR2093959
SRR2093960
SRR2093961
SRR2089754

建立一个sra文件夹，使用prefetch将数据下载到该文件夹下

1 2	mkdir sra cat SRR_Acc_List.txt \| while read id; do prefetch -O sra $id ; done

然后将sra文件转成fastq格式

1 2	mkdir -p 00-raw-data cat SRR_Acc_List.txt \| while read id; do fasterq-dump -p -3 -O 00-raw-data sra/$id.sra ; done &

原始数据质控

原文只说了High quality (QV > 25) 的read用于比对，并没有提到具体的质控的软件。我们这里用fastp做质控，因为fastp跑完之后也会输出一个网页，所以我也就懒得用FastQC再跑一次

1
2
3

mkdir -p 01-clean-data
mkdir -p qc
cat SRR_Acc_List.txt | while read id; do fastp -w 4 -i 00-raw-data/$id.sra_1.fastq -I 00-raw-data/$id.sra_2.fastq -o 01-clean-data/$id.sra_1.fastq.gz -O 01-clean-data/$id.sra_2.fastq.gz --html qc/$id.html --json qc/$id.json ; done &

重申一个观点，对于普通转录组，由于测序质量过关，比对软件也能处理开头的接头序列，我基本不做质控。

建立索引

文章说它用到水稻参考基因组和注释文件下载自ftp://ftp.plantbiology.msu.edu/pub/data/EukaryoticProjects/osativa/annotationdbs/pseudomolecules/version_7.0/all.dir/, 而有趣的是这个链接我打不开，目前可用的下载地址其实是: http://rice.plantbiology.msu.edu/pub/data/Eukaryotic_Projects/o_sativa/annotation_dbs/pseudomolecules/version_7.0/all.dir/

mkdir reference && cd reference
# 参考基因组
wget http://rice.plantbiology.msu.edu/pub/data/Eukaryotic_Projects/o_sativa/annotation_dbs/pseudomolecules/version_7.0/all.dir/all.con
# 注释文件
wget http://rice.plantbiology.msu.edu/pub/data/Eukaryotic_Projects/o_sativa/annotation_dbs/pseudomolecules/version_7.0/all.dir/all.gff3

此外文章用的是TopHat v1.3.3和Bowtie v0.12.9 默认参数做的比对，我对此不作任何评价。我们这里用HISAT2建立索引

1	hhisat2-build -p 10 reference/all.con reference/Osativa

序列比对

用HISAT2将过滤后的fastq文件和参考基因组进行比对，得到下游分析所需要的BAM文件

1
2

mkdir -p 02-read-align
cat SRR_Acc_List.txt | while read id; do hisat2 -p 10 -x reference/Osativa -1 01-clean-data/$id.sra_1.fastq.gz -2 01-clean-data/$id.sra_2.fastq.gz --new-summary --summary-file qc/$id.ht2.txt | samtools sort -@ 4 > 02-read-align/$id.bam ; done &

表达量定量

文章使用的Cufflinks v1.3.0对基因的表达量进行定量，使用FPKM进行了标准化。我这里也不做评价，我们用featureCounts进行表达量定量

先将GFF3转成GTF格式

1	gffread reference/all.gff3 -T -o reference/Osativa.gtf

然后用featureCoutns

1
2
3

mkdir -p 03-gene-count
featureCounts -a reference/Osativa.gtf -o 03-gene-count/gene_counts 02-read-align/*.bam -T 10 -g gene_id
mv 03-gene-count/gene_counts.summary qc/gene_counts.summary

最终的结果，就可以读取到R语言里面进行分析了

MultiQC

上面每一步的中间文件都被我放到了qc文件夹里，因此可以用multiqc进行整合

差异分析

由于前面几步的运算时间比较久，我这里直接准备了表达量矩阵，下载方式为百度网盘链接：https://pan.baidu.com/s/1c_mdGlDggV2DH79H-X7J_g 提取码：jmi5

差异表达分析用但是精确检验(exact test)。

之前的SRR编号在分析的时候最好转成实际的样本编号，方便分析的时候查看

编号	品种	处理
SRR2089751.sra	Sahabhagidhan	0h control（易感品种对照）
SRR2089753.sra	Sahabhagidhan	6h冷害处理
SRR2089754.sra	Sahabhagidhan	12h冷害处理
SRR2089756.sra	Sahabhagidhan	24h冷害处理
SRR2093937.sra	Sahabhagidhan	48h冷害处理
SRR2093938.sra	Sahabhagidhan	48h冷害处理后又24h恢复
SRR2093939.sra	Geetanjali	0h conrol
SRR2093948.sra	Geetanjali	6h冷害处理
SRR2093955.sra	Geetanjali	12h冷害处理
SRR2093959.sra	Geetanjali	24h冷害处理
SRR2093960.sra	Geetanjali	48h冷害处理
SRR2093961.sra	Geetanjali	48h冷害处理后又24h恢复

Sahabhagidhan是不耐冷的品种(cold susceptible variety, CSV) ,而Geetanjali则是耐冷品种 (CTV, cold tolerant variety) 。

我们这里分析的时候把两个物种在6, 12, 24, 48时间点的样本作为重复，仅仅比较不同物种在冷害处理后的差异基因，而不分析冷害处理后随时间点变化的基因。

加载数据:

df <- read.table("gene_counts", header = T,
                 sep = "\t",
                 stringsAsFactors = F)
expr_df <- df[,c(1,8,9,10,11,14,15,16,17)]
colnames(expr_df) <- c("geneid",
                       "CSV_1","CSV_2","CSV_3","CSV_4",
                       "CTV_1","CTV_2","CTV_3","CTV_4")

DESeq2差异分析

library(DESeq2)
coldata <- data.frame(condition=factor(c("CSV","CSV","CSV","CSV",
                      "CTV","CTV","CTV","CTV")))
dds <- DESeqDataSetFromMatrix(expr_df,
                       colData = coldata,
                       design = ~ condition,
                       tidy = TRUE)
dds <- DESeq(dds)

PCA比较下组内差异是不是比组间差异小

1 2	rld <- rlog(dds) plotPCA(rld)

PCA

DESeq2的差异表达分析，只需要一行代码

1	res <- results(dds)

用ggplot2画个火山图,

plot_df <- data.frame(res)
plot_df <- cbind(geneid=row.names(plot_df),
                        plot_df)

plot_df <- plot_df[!is.na(plot_df$padj), ]

ggplot(data=plot_df, aes(x=log2FoldChange, 
                      y =-log10(padj))) +
  geom_point(alpha=0.8, size=1.2)+
  labs(title="Volcanoplot", x="log2 (fold change)",y="-log10 (q-value)")+
  theme(plot.title = element_text(hjust = 0.4))+
  geom_hline(yintercept = -log10(0.05),lty=4,color="red")+
  geom_hline(yintercept = -log10(0.01),lty=4,color="blue")+
  geom_vline(xintercept = c(1,-1),lty=4,alpha=0.8,color="blue")+
  geom_vline(xintercept = c(2,-2),lty=4,alpha=0.8,color="red")+
  theme_bw()+
  theme(panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(colour = "black"))

火山图

火山图的作用，一是检查你的差异分析是否合理，也就是大部分基因都是没有变化，只有一部分是差异基因，第二个是可以帮你确定阈值。

我发现即便是LogFoldChange > 2（相差4倍）, p.adj < 0.01，依旧是有很多差异基因。所以我这里以这个为标准进行筛选，

res_sub <- subset(res, abs(log2FoldChange) > 2 & padj <0.01)
gene_up <- row.names(res_sub[res_sub$log2FoldChange > 0, ])
length(gene_up) # 1024
gene_down <- row.names(res_sub[res_sub$log2FoldChange < 0, ])
length(gene_down) # 702

其实这个时候你会有一个问题，所谓的基因上调和下调是谁相对于谁呢？很简单，直接把差异的基因的标准化矩阵提取出来即可。

1
2
3

gene_counts <- counts(dds, normalized = TRUE)
gene_counts <- gene_counts[row.names(res_sub), ]
head(gene_counts)

当你看到具体的数值时

                    CSV_1     CSV_2      CSV_3      CSV_4     CTV_1
LOC_Os01g01430   1.739379   0.00000   5.971535   8.188634  31.70621
                   CTV_2     CTV_3    CTV_4
LOC_Os01g01430  29.60817  17.19996  10.9967

你自然就知道了CTV vs CSV里，CSV是作为分母，而CTV是分子。

log2 fold change (MLE): condition CTV vs CSV 
Wald test p-value: condition CTV vs CSV 
DataFrame with 6 rows and 6 columns
                       baseMean    log2FoldChange             lfcSE
                      <numeric>         <numeric>         <numeric>
LOC_Os01g01430 13.1763234817146  2.57839743463268 0.925763093663462

使用clusterProfiler做GO/KEGG富集分析

文章用agriGo对DEG进行分类，用singular enrichment analysis (SEA) 做GO富集分析，进一步的注释用Web Gene Ontology Annotation Plot (WEGO) 进行画图。做KEGG富集分析时，用的是KEGG GENES database KAAS (KEGG Automatic Annotation Server) 对差异基因进行功能注释.

这里我也不做评价，我们此处用Y叔的clusterProfiler做GO/KEGG富集分析。

在做GO富集分析的时候，我遇到了一个难题，我发现Bioconductor上并没有水稻的物种注释包。解决这个问题有两种方法，

一种是利用AnnotationHub在线检索抓取OrgDb, 但是这些包都是用ENTREZID，你需要先将RAP-DB或者MSU转为ENTREZID才行
另一种则是我构建好一个物种注释包，你们只要安装就可以继续使用enrichGO，参考https://github.com/xuzhougeng/org.Osativa.eg.db

参考我的GitHub页面安装，https://github.com/xuzhougeng/org.Osativa.eg.db

1
2
3

install.packages("https://github.com/xuzhougeng/org.Osativa.eg.db/releases/download/v0.01/org.Osativa.eg.db.tar.gz", 
    repos = NULL, 
    type="source")

目前这个物种包比较简单，就只能做GO富集分析，以及MSU和RAP-DB之间的ID转换

然后进行GO富集分析

library(clusterProfiler)
library(org.Osativa.eg.db)

org <- org.Osativa.eg.db
ego_up <- enrichGO(gene_up,
                OrgDb = org,
                keyType = "GID",
                pAdjustMethod = "none",
                ont="BP")
p1 <- dotplot(ego_up)

ego_down <- enrichGO(gene_down,
                OrgDb = org,
                keyType = "GID",
                pAdjustMethod = "none",
                ont="BP")
p2 <- dotplot(ego_down)

library(cowplot)
plot_grid(p1,p2)

dotplot

除了点图，我们还可以用clusterProfiler画下面这个图，展现每个基因在不同GO通路的倍数变化

fc <- plot_df$log2FoldChange
names(fc) <-  plot_df$geneid

heatplot(ego_down, foldChange = fc)

heatplot

讲真，GO富集分析给的信息真的是少，我以为是我的分析出现了问题，所以我去看了文章里GO分析中BP部分，发现信息是一样的少。

GO富集分析做完之后，KEGG富集分析怎么做呢？参考我之前的这篇https://mp.weixin.qq.com/s/UnUPVoaMpfJWCQEpkmdTWA

因为KEGG要求的输入编号是RAP-DB，所以我们需要进行一次ID转换，在编号中加上”-01”，将编号中的”g”改成”t”, 就可以做KEGG富集分析了

rap_id <- mapIds(x = org, 
                 keys = gene_up, 
                 column = "RAP",
                 "GID")
rap_id <- paste0(rap_id[!is.na(rap_id)], "-01")
rap_id <- gsub("g","t",rap_id)
ekegg <- enrichKEGG(rap_id, organism = "dosa", pAdjustMethod = "none")

没成想，KEGG的分析结果中只得到了一个结果, 所以也懒得放气泡图了。也难怪文章只是将差异基因按照KEGG进行了分类，而没有做富集。我这里也用Y叔的包对上调差异基因按照KEGG归类

rice_kegg <- clusterProfiler::download_KEGG("dosa")

kegg_df <- rice_kegg$KEGGPATHID2EXTID
kegg_df <- kegg_df[kegg_df$to %in% rap_id,]
kegg_df <- merge(kegg_df, rice_kegg$KEGGPATHID2NAME,
                 by.x="from",by.y="from")

kegg_class <- as.data.frame(sort(table(kegg_df$to.y), decreasing = T)[1:10])
colnames(kegg_class) <- c("pathway","times")
ggplot(kegg_class,aes(x=pathway, y = times)) +
  geom_bar(fill="#ca0020",stat="identity") + coord_flip() +
  theme_bw() + geom_text(aes(y = times+1, label = times))

KEGG分类1

当然上下调一起展示也不难，就是代码会长一点

rice_kegg <- clusterProfiler::download_KEGG("dosa")
# 获取上调的KEGG PATHWAY
up_rap_id <- mapIds(x = org, 
                 keys = gene_up, 
                 column = "RAP",
                 "GID")
up_rap_id <- paste0(up_rap_id[!is.na(up_rap_id)], "-01")

up_rap_id <- gsub("g","t",up_rap_id)

up_df <- rice_kegg$KEGGPATHID2EXTID
up_df <- kegg_df[kegg_df$to %in% up_rap_id,]
up_df <- merge(kegg_df, rice_kegg$KEGGPATHID2NAME,
                 by.x="from",by.y="from")

# 获取下调的KEGG PATHWAY
down_rap_id <- mapIds(x = org, 
                    keys = gene_down, 
                    column = "RAP",
                    "GID")
down_rap_id <- paste0(down_rap_id[!is.na(down_rap_id)], "-01")
down_rap_id <- gsub("g","t",down_rap_id)
down_df <- rice_kegg$KEGGPATHID2EXTID
down_df <- down_df[down_df$to %in% down_rap_id,]
down_df <- merge(down_df, rice_kegg$KEGGPATHID2NAME,
               by.x="from",by.y="from")

# 统计
kegg_class_up <- as.data.frame(sort(table(up_df$to.y), 
                                 decreasing = T)[1:10])
kegg_class_down <- as.data.frame(sort(table(down_df$to.y), 
                                    decreasing = T)[1:10])
#合并
kegg_class <- rbind(kegg_class_up, kegg_class_down)
colnames(kegg_class) <- c("pathway","times")
kegg_class$source <- rep(c("up","down"),
                         times=c(nrow(kegg_class_up),nrow(kegg_class_down)))
#画图
ggplot(kegg_class,aes(x=pathway, y = times)) +
  geom_bar(aes(fill=source),stat="identity",position = "dodge") + 
  scale_fill_manual(values = c(up="#ca0020",down="#2b83ba")) +
  coord_flip() +
  theme_bw()

KEGG分类2

这里这展示了代码，具体如何讲故事，如何调整参数，其实就靠作者自由发挥了。

参考资料

水稻两大权威注释组织

RAP-DB: https://rapdb.dna.affrc.go.jp/
RGAP: http://rice.plantbiology.msu.edu/
clusterPrifiler文档: https://bioconductor.org/packages/release/bioc/vignettes/clusterProfiler/inst/doc/clusterProfiler.html
DESeq2文档: http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html

2019-08-22

其他

如何让你的GitHub能被引用

你需要先登陆到https://zenodo.org/，然后需要创建一个账号，这个账号一定要用GitHub账号进行授权。之后你发现一个Using GitHub，点击Check out

Step1

然后你就会发现自己的repository都被列出来了

Step2

之后将你想要被引用的仓库的按钮从off变成ON

Step3

接着到你想被引用的仓库里, 我们要创建要给release，这样子就有了版本号，如此别人才能知道自己用的版本是哪一版。

Step4

我们就可以创建release，可以设置合适的版本好，比如说我这边设置为0.01

Step5

我们就拥有了一个释放版本

Step6

之后回到Zenodo，点击Upload

Step7

就会出现你刚才释放的软件

Step8

点击它，你就拥有了可以引用的DOI号

Step9

参考资料：https://guides.github.com/activities/citable-code/

2019-08-21

生信软件工具箱

MECAT2:三代高效组装工具

MECAT2: 三代高效组装工具

MECAT2是三代测序数据PacBio的高效组装工具，是之前MECAT的改进版, 修复了之前的很多bug, 使用基于string graph的fsa替换了之前的mecat2canu。

MECAT2由4个模块组成:

mecat2pw: SMART reads快速和准确地配对比对工具
mecat2ref: SMART reads的参考基因组比对工具
mecat2cns: 基于配对的重叠对存在噪音的reads进行纠错
fsa: 基于string graph的组装工具

目前MECAT2只支持PacBio数据，对于Nanopore数据，肖老师他们开发了NECAT

安装

MECAT2软件安装非常简单，不存在上一代MECAT的安装难问题了

mkdir -p ~/opt/biosoft
cd ~/opt/biosoft
git clone https://github.com/xiaochuanle/MECAT2.git
cd MECAT2
make -j 4

唯一要注意的是尽量要用新版的perl，经测试，Perl 5.26 是可以正常运行脚本。

最后加入到环境变量PATH中

1	export PATH=~/opt/biosoft/MECAT2/Linux-amd64/bin:$PATH

软件使用

我们以 E. coli 数据集作为案例，介绍MECAT2应该如何使用。

第一步，我们要下载数据

mkdir -p ~/ecoli
cd ~/ecoli
wget http://gembox.cbcb.umd.edu/mhap/raw/ecoli_filtered.fastq.gz
gunzip ecoli_filtered.fastq.gz

注意: 目前MECAT2还不支持gz压缩文件，如果你直接用gz作为输入，会提示core dumped。

MECAT2不是根据文件名来确定输入的序列格式, 也就是如果你把 ecoli_filtered.fastq 命名为 ecoli_filtered.fasta, 组装也不会出错。

第二步, 准备模板的config文件

1	mecat.pl config ecoli_config_file.txt

使用vim修改config文件

PROJECT=ecoli
RAWREADS=ecoli_filtered.fastq
GENOME_SIZE=4800000
THREADS=4
MIN_READ_LENGTH=2000
CNS_OVLP_OPTIONS=""
CNS_OPTIONS="-r 0.6 -a 1000 -c 4 -l 2000"
TRIM_OVLP_OPTIONS="-B"
ASM_OVLP_OPTIONS="-n 100 -z 10 -b 2000 -e 0.5 -j 1 -u 0 -a 400"
FSA_OL_FILTER_OPTIONS="--max_overhang=-1 --min_identity=-1"
FSA_ASSEMBLE_OPTIONS=""
USE_GRID=false
CLEANUP=0
CNS_OUTPUT_COVERAGE=30
GRID_NODE=0

第三步: 原始数据纠错

1	mecat.pl correct ecoli_config_file.txt

第四步: 对纠错后的reads进行组装

1	mecat.pl assemble ecoli_config_file.txt

输出结果:

纠错后的reads： ecoli/1-consensus/cns_reads.fasta.
最长30X纠错后用于trimming的reads: ecoli/1-consensus/cns_final.fasta.
trimmed reads: ecoli/2-trim_bases/trimReads.fasta
组装的contigs: ecoli/4-fsa/contigs.fasta

最终的组装结果为

Count: 1
Tatal: 4630437
Max: 4630437
Min: 4630437
N25: 4630437
L25: 1

能够完美的就装出一条序列

配置文件介绍

MECAT2的组装过程的参数都在配置文件中, 也就是之前的ecoli_config_file.txt。我按照不同部分对这些参数进行介绍

基本参数

这些参数属于改起来不怎么纠结的参数

PROJECT=ecoli: 项目名，之后的所有输出文件都在该项目名的目录下
RAWREADS=: 原始数据的位置，可以是Fasta或者是Fastq， H5格式要先转成Fasta (参考 Input Format).
GENOME_SIZE=: 基因组大小，单位是bp，如果是120M，那么就是120000000.
THREADS=: 线程数
MIN_READ_LENGTH=: 用于纠错和trim的reads的最低长度. 数据质量好，就长一点，比如说1000到2000
USE_GRID=false: 是否有多个计算节点
CLEANUP=0: 运行结束后删除MECAT2的中间文件, 大基因组的临时文件很大，所以要设置为1.
GRID_NODE=0, 当USE_GRID=1时，设置用到的计算节点数，如果是单节点服务器，不需要设置。

纠错和修剪阶段

纠错阶段(correct stage)和修剪阶段(trimming stage)，MECAT2调用的mecat2pw和mecat2cns，与之相关的配置如下

CNS_OVLP_OPTIONS="", 在纠错阶段是检测候选重叠的参数, 会传给mecat2pw
CNS_OPTIONS=""：原始reads纠错参数，会传递给mecat2cns,
TRIM_OVLP_OPTIONS="": 在trim阶段，用于检测重叠的参数,会传给v2asmpm

但是我发现v2asmpm没有文档说明，和软件开发者讨论之后，给出的说明如下

V2asmpm -Pworkpath -Tt -Sx -Ey -B

-Pworkpath 工作目录是workpath
-Tt 用t个CPU线程
-Sx -Ey 这两个参数一起表示只计算第x到第y卷(包括第y卷), 通常工作目录下会有000001.fasta, 000002.fasta这样的分卷
-B 如果有这个选项, 就输出二进制格式的比对结果; 如果没有这个选项, 就输出文本格式的比对结果

因此用默认的-B就行了。

组装阶段

组装阶段相关的参数如下

CNS_OUTPUT_COVERAGE=30: 选择多少覆盖度的最长纠错后reads进行trim和组装。举个例子，4.8M的 E. coli 30X, 是144MB。一般选择20~40之间就行，会传递给mecat2elr
ASM_OVLP_OPTIONS="": 在组装阶段，用于检测重叠的参数，传给v2asmpm.sh
FSA_OL_FILTER_OPTIONS="", 过滤重叠的参数，传递给fsa_ol_filter
FSA_ASSEMBLE_OPTIONS="", 组装trimm reads的参数, 传给v2asmpm

这些参数调整起来都很复杂，只能建议看函数帮助文档了。举个例子FSA_OL_FILTER_OPTIONS="--max_overhang=-1 --min_identity=-1"根据阅读fsa_ol_filter的帮助发现，-1表示让程序自己决定

而ASM_OVLP_OPTIONS="-n 100 -z 10 -b 2000 -e 0.5 -j 1 -u 0 -a 400"，根据和软件开发者的讨论，是不需要额外调整，用默认的即可。

MECAT2是基于比对纠错的组装工具，组装速度上比同类型的Canu和Falcon都快。仅看组装结果的N50参数，MECAT2结果和Canu, Falcon结果都是差不多的。

考虑到BUSCO完整度，Canu是三者最高。不过用三代数据和二代数据进行几轮Polish后会提高到Canu相同水平。三个软件上的BUSCO（embryophyta_odb10）值如下:

1
2
3

Canu: C:98.8%[S:96.3%,D:2.5%],F:0.5%,M:0.7%,n:1375
MECAT2: C:93.5%[S:91.6%,D:1.9%],F:4.6%,M:1.9%,n:1375
Falcon: C:96.8%[S:95.1%,D:1.7%],F:2.0%,M:1.2%,n:1375

因此，如果你原本你是用Falcon做组装，那么可以改成MECAT2了，效果差不多，速度还更快了。对于越大的物种，在有限的资源下，更推荐用MECAT2。

参考资料

https://github.com/xiaochuanle/MECAT2

2019-08-20

生信软件工具箱

「三代组装」使用新版Falcon进行三代测序基因组组装

这里的新版指的是PacBio公司在2018年9月发布pb-assembly, 而这篇文章是在2018年9月30日发的。

今年早些时候在参加三代培训时，听说PacBio会在今年对Falcon进行一些改变。前几天我在读 readthedocs上的Falcon文档时，发现了文档页面上方出现了这样两栏提醒

其中pb_assembly就是新的FALCON组装套装在GitHub上的使用文档，经过了几天的探索，我对它有了一点了解，写一篇教程作为官方文档的一些补充吧。

新版亮点

整合了串联重复序列和遍在重复序列的屏蔽（之前没有这一步）
以GFA格式存放graph文件，后续可以用Bandage进行可视化
通过算法和性能优化，提高了Associate Contigs的准确性
分析流程的性能优化

软件安装和数据准备

Falcon终于拥抱了bioconda, 这也就意味着我们再也不需要用到他们原本笨拙的安装脚本，浪费时间在安装软件上。

conda create -n pb-assembly  pb-assembly
source activate pb-assembly
# 或者
conda create -p ~/opt/biosoft/pb-assembly pb-assembly
source activate ~/opt/biosoft/pb-assembly

这里使用https://pb-falcon.readthedocs.io/en/latest/tutorial.html上的所用的E. coli数据集

1
2

wget https://downloads.pacbcloud.com/public/data/git-sym/ecoli.m140913_050931_42139_c100713652400000001823152404301535_s1_p0.subreads.tar.gz
tar -xvzf ecoli.m140913_050931_42139_c100713652400000001823152404301535_s1_p0.subreads.tar.gz

解压缩后的文件夹里有三个300M的Fasta文件, 将他们的实际路径记录到input.fofn中

1
2
3

ecoli.1.fasta
ecoli.2.fasta
ecoli.3.fasta

准备配置文件

为了进行组装，需要准备一个配置文件。我的配置文件为fc_run.cfg，内容如下。你们可以先预览一下，后面看我的解释说明。

#### Input
[General]
input_fofn=input.fofn
input_type=raw
pa_DBdust_option=
pa_fasta_filter_option=pass
target=assembly
skip_checks=False
LA4Falcon_preload=false

#### Data Partitioning
pa_DBsplit_option=-x500 -s50
ovlp_DBsplit_option=-x500 -s50

#### Repeat Masking
pa_HPCTANmask_option=
pa_REPmask_code=1,100;2,80;3,60

####Pre-assembly
genome_size=0
seed_coverage=20
length_cutoff=3000    
pa_HPCdaligner_option=-v -B128 -M24
pa_daligner_option=-e.7 -l1000 -k18 -h80 -w8 -s100
falcon_sense_option=--output-multi --min-idt 0.70 --min-cov 2 --max-n-read 800
falcon_sense_greedy=False

####Pread overlapping
ovlp_daligner_option=-e.96 -l2000 -k24 -h1024 -w6 -s100
ovlp_HPCdaligner_option=-v -B128 -M24

####Final Assembly
overlap_filtering_setting=--max-diff 100 --max-cov 300 --min-cov 2
fc_ovlp_to_graph_option=
length_cutoff_pr=2000

[job.defaults]
job_type=local
pwatcher_type=blocking
JOB_QUEUE=default
MB=32768
NPROC=6
njobs=32
submit = /bin/bash -c "${JOB_SCRIPT}" > "${JOB_STDOUT}" 2> "${JOB_STDERR}"

[job.step.da]
NPROC=4
MB=32768
njobs=32
[job.step.la]
NPROC=4
MB=32768
njobs=32
[job.step.cns]
NPROC=8
MB=65536
njobs=5
[job.step.pla]
NPROC=4
MB=32768
njobs=4
[job.step.asm]
NPROC=24
MB=196608
njobs=1

根据注释信息，文件分为”input”, “Data partitioning”, “Repeat Masking”, “Pre-assembly”, “Pread overlapping”, “Final Assembly”, 以及最后的任务调度部分，让我们分别看下这里面的内容

输入(Input)

#### Input
[General]
input_fofn=input.fofn
input_type=raw
pa_DBdust_option=
pa_fasta_filter_option=pass
target=assembly
skip_checks=False
LA4Falcon_preload=false

输入这里参数比较简单，基本不需要做任何改动，除了 pa_fasta_filter_options，用于处理一个ZMW（测序翻译孔）有多条subread时，到底选择哪一条的问题。

“pass”: 不做过滤，全部要。
“streamed-median”: 表示选择大于中位数的subread
“streamed-internal-median”: 当一个ZMW里的subread低于3条时选择最长，多于单条则选择大于中位数的subread

0.01版本pb-assembly的pa_DBdust_option有一个bug，也就是里面的参数不会传递给DBdust, DBdust是对read进行soft-masking，一般都用默认参数，因此这个bug问题不大。

数据分配(Data Partitioning)

1 2	pa_DBsplit_option=-x500 -s50 ovlp_DBsplit_option=-x500 -s50

这部分的设置会将参数传递给DBsplit，将数据进行拆分多个block，后续的并行计算都基于block。-s 50表示每个block大小为50M。这适用于基因组比较小的物种，如果是大基因组则应该设置为-s 200或者-s 400

重复屏蔽(Repeat Masking)

1
2
3

#### Repeat Masking
pa_HPCTANmask_option=
pa_REPmask_code=1,100;2,80;3,60

屏蔽重复序列可以在不损失组装准确性的同时，提高后续组装的overlap/daligner步骤10~20倍速度，见Detecting and Masking Repeats.

pa_HPCTANmask_option的参数会传给串联重复步骤的HPCTANmask, 而pa_REPmask_code很复杂，它分为三次迭代，因此这里1:100;2,80;3,60 就表示第一次迭代检测每个block中出现超过100次的序列，第二次迭代将2个block合并一起检测超过80次的序列，第三次将3个block进行合并检测超过60次的序列。

https://github.com/PacificBiosciences/pbbioconda/issues/20

预组装（纠错）pre-assembly

####Pre-assembly
genome_size=0
seed_coverage=20
length_cutoff=3000    
pa_HPCdaligner_option=-v -B128 -M24
pa_daligner_option=-e.7 -l1000 -k18 -h80 -w8 -s100
falcon_sense_option=--output-multi --min-idt 0.70 --min-cov 2 --max-n-read 800
falcon_sense_greedy=False

当length_cutoff=-1时，设置genome_size和seed_coverage会自动计算要过滤的序列。否则是过滤低于一定长度的read。

pa_HPCdalinger_option参数不需要调整，-M表示每个进程的内存为24G，一般200M的block对应16G。

pa_daligner_option的参数比较重要：

-e:错误率，低质量序列设置为0.70，高质量设置为0.80。值越高避免单倍型的坍缩
-l: 最低overlap的长度，文库比较短时为1000，文库比较长为5000.
-k: 低质量数据为14，高质量数据为18
-h: 表示完全match的k-mer所覆盖的碱基数。和-l, -e有关，越大越严格。预组装时最大也不要超过最低overlap长度的1/4. 最低就设置为80

纠错后相互比对

1
2
3

####Pread overlapping
ovlp_daligner_option=-e.96 -l2000 -k24 -h1024 -w6 -s100
ovlp_HPCdaligner_option=-v -B128 -M24

和上面的参数类似，但是-e的范围调整为0.93-0.96，-l范围调整为1800-6000， -k调整为 18-24

最后组装

1
2
3

overlap_filtering_setting=--max-diff 100 --max-cov 300 --min-cov 2
fc_ovlp_to_graph_option=
length_cutoff_pr=2000

这里的参数就可以随便调整了，因为这一步速度很快。例如length_cutoff_pr就可以从2000，提高到15000.

最后还有一部分是任务投递系统，如果是单节点运行，需要注意设置 njobs，这是同时投递的任务数。假如你将[job.step.cns]按照如下的方式设置，那么同时会出现 $ 8 \times 50 = 400 $ 个任务，如果你的内存只有128G，运行一段时间后你的所有内存就会被耗尽，那么基本上你就只能重启服务器了。

[job.step.cns]
NPROC=8
MB=65536
njobs=50

运行结果

用上述的配置文件，以fc_run fc_run.cfg运行后，最后的2-asm-falcon/p_ctg.fa的序列数有4条，最长为4,685,024, 之后我将length_cutoff_pr调整为15k，2-asm-falcon/p_ctg.fa序列只有一条，长度为4,638,411

下载asm.gfa到本地，用Bandage可视化，可以发现组装效果不错。

Bandage

2019-08-18

生信软件工具箱

Cicero:单细胞共开放分析

Cicero是一个单细胞染色质可及性实验可视化R包。Cicero的主要功能就是使用单细胞染色质可及性数据通过分析共开放去预测基因组上顺式调节作用(cis-regulatory interactions)，例如增强子和启动子。Cicero能够利用染色质开放性进行单细胞聚类，排序和差异可及性分析。关于Cicero的算法原理，参考他们发表的文章.

简介

Cicero的主要目标是使用单细胞染色质可及性数据去预测基因组上那些在细胞核中存在物理邻近的区域。这能够用于鉴定潜在的增强子-启动子对，对基因组顺式作用有一个总体了解。

由于单细胞数据的稀疏性，细胞必须根据相似度进行聚合，才能够对数据中的多种技术因素进行可靠纠正。

最终，Cicero给出了”Cicero 共开放”得分，分数在-1和1之间，用于评估用于给定距离中每个可及peak对之间的共开放性，分数越高，越可能是共开放。

此外，Monocle3的框架让Cicero能对单细胞ATAC-seq实验完成Monocle3上的聚类，拟时间等分析。

Cicero主要提供了两种核心分析:

构建和分析顺式调控网络.
常规单细胞染色质可及性分析:

Cicero安装和加载

Cicero运行在R语言分析环境中。尽管能够从Bioconductor安装，但是推荐从GitHub上安装。

install.packages("BiocManager")
BiocManager::install(c("Gviz", "GenomicRanges", "rtracklayer"))
install.packages("devtools")
devtools::install_github("cole-trapnell-lab/cicero-release", ref = "monocle3")

加载Cicero

1	library(cicero)

数据加载

Cicero以cell_data_set(CDS)对象存放数据，该对象继承自Bioconductor的SingleCellExperiment对象。我们可以用下面三个函数来操纵数据

fData: 获取feature的元信息, 这里的featureI指的是peak
pData: 获取cell/sample的元信息
assay: 获取cell-by-peak的count矩阵

除了CDS外，还需要基因坐标信息和基因注释信息。

关于坐标信息，以mouse.mm9.genome为例，它是一个数据框，只有两列，一列是染色体编号，另一列对应的长度

1	data("mouse.mm9.genome")

基因组注释信息可以用rtracklayer::readGFF加载

temp <- tempfile()
download.file("ftp://ftp.ensembl.org/pub/release-65/gtf/mus_musculus/Mus_musculus.NCBIM37.65.gtf.gz", temp)
gene_anno <- rtracklayer::readGFF(temp)
unlink(temp)

Bioconductor也有现成的TxDb，可以从中提取注释信息。

简单稀疏矩阵格式

以一个测试数据集介绍如何加载简单稀疏矩阵格式

1 2	temp <- textConnection(readLines(gzcon(url("http://staff.washington.edu/hpliner/data/kidney_data.txt.gz")))) cicero_data <- read.table(temp)

这里的简单稀疏矩阵格式指的是数据以三列进行存放，第一列是peak位置信息，第二列是样本信息，第三列则是样本中有多少fragment和该peak重叠。

1	chr1_4756552_4757256 GAATTCGTACCAGGCGCATTGGTAGTCGGGCTCTGA 1

对于这种格式，Cicero提供了一个函数make_atac_cds, 用于构建一个有效的cell_data_set对象，用于下游分析，输入既可以是一个数据框，也可以是一个文件路径。

1	input_cds <- make_atac_cds(cicero_data, binarize = TRUE)

10X scATAC-seq 数据

如果数据已经经过10X scATAC-seq处理，那么结果中的filtered_peak_bc_matrix里就有我们需要的信息

加载cell-by-peak count矩阵，

# read in matrix data using the Matrix package
indata <- Matrix::readMM("filtered_peak_bc_matrix/matrix.mtx") 
# binarize the matrix
indata@x[indata@x > 0] <- 1

加载cell元数据

# format cell info
cellinfo <- read.table("filtered_peak_bc_matrix/barcodes.tsv")
row.names(cellinfo) <- cellinfo$V1
names(cellinfo) <- "cells"

加载peak元数据

# format peak info
peakinfo <- read.table("filtered_peak_bc_matrix/peaks.bed")
names(peakinfo) <- c("chr", "bp1", "bp2")
peakinfo$site_name <- paste(peakinfo$chr, peakinfo$bp1, peakinfo$bp2, sep="_")
row.names(peakinfo) <- peakinfo$site_name

为cell-by-peak的count矩阵增加行名(peak元数据)和列名(cell元数据)

1 2	row.names(indata) <- row.names(peakinfo) colnames(indata) <- row.names(cellinfo)

然后用new_cell_data_set根据peak元数据，cell元数据构建出一个cell_data_set对象

# make CDS
input_cds <-  suppressWarnings(new_cell_data_set(indata,
    cell_metadata = cellinfo,
    gene_metadata = peakinfo))

之后用detect_genes统计，对于每个基因在多少细胞中的表达量超过了给定阈值。

1	input_cds <- monocle3::detect_genes(input_cds)

过滤没有细胞开放的peak

1	input_cds <- input_cds[Matrix::rowSums(exprs(input_cds)) != 0,]

构建顺式作用网络

运行Cicero

构建Cicero CDS对象

因为单细胞染色质开放数据特别稀疏，为了能够比较准确的估计共开放得分，需要将一些比较相近的细胞进行聚合，得到一个相对比较致密的count数据。Cicero使用KNN算法构建细胞的相互重叠集。而细胞之间距离关系则是根据降维后坐标信息。降维方法有很多种，这里以UMAP为例

input_cds <- detect_genes(input_cds)
input_cds <- estimate_size_factors(input_cds)
# PCA或LSI线性降维
input_cds <- preprocess_cds(input_cds, method = "LSI")
# UMAP非线性降维
input_cds <- reduce_dimension(input_cds, reduction_method = 'UMAP', 
                              preprocess_method = "LSI")

此时用reducedDimNames(input_cds)就会发现它多了LSI和UMAP两个信息，我们可以用reducedDims来提取UMAP坐标信息。

1	umap_coords <- reducedDims(input_cds)$UMAP

make_cicero_cds函数就需要其中UMAP坐标信息。

注1: 假如你已经有了UMAP的坐标信息，那么就不需要运行preprocess_cds和reduce_dimension，直接提供UMAP坐标信息给make_cicero_cds就可以了。

注2: umap_coords的行名需要和pData得到表格中的行名一样, 即all(row.names(umap_coords) == row.names(pData(input_cds)))结果为TRUE。

1	cicero_cds <- make_cicero_cds(input_cds, reduced_coordinates = umap_coords)

此时的cicero_cds是没有UMAP降维信息。

运行Cicero

Cicero包的主要功能就是预测基因组上位点之间的共开放. 有两种方式可以获取该信息

run_cicero: 一步到位。默认参数适用于人类和小鼠，但是其他物种，参考此处
单独运行每个函数，这样子获取中间的信息。

#测试部分数据
data("mouse.mm9.genome")
sample_genome <- subset(mouse.mm9.genome, V1 == "chr2")
sample_genome$V2[1] <- 10000000
conns <- run_cicero(cicero_cds, sample_genome, sample_num = 2) 
nrow(conns) # 212302
## 全部数据, 时间真久
conns <- run_cicero(cicero_cds, mouse.mm9.genome, sample_num = 2)
nrows(conns) #59741738

其中conns存放的是peak之间共开放得分。

Cicero Connection可视化

有了peak之间共开放的可能性得分后，就可以画一些很酷的图了。用到的函数是plot_connections, 这个函数有很多参数，所以最终呈现结果取决于你的参数调整。

获取基因组注释信息

temp <- tempfile()
download.file("ftp://ftp.ensembl.org/pub/release-65/gtf/mus_musculus/Mus_musculus.NCBIM37.65.gtf.gz", temp)
gene_anno <- rtracklayer::readGFF(temp)
unlink(temp)

gene_anno$chromosome <- paste0("chr", gene_anno$seqid)
gene_anno$gene <- gene_anno$gene_id
gene_anno$transcript <- gene_anno$transcript_id
gene_anno$symbol <- gene_anno$gene_name

除了从GFF文件中提取外，还可以利用txdb提取转录本位置信息

举例子，下面的函数就是绘制染色体上"chr2:9773451-9848598的共开放情况，

plot_connections(conns, "chr2", 9773451, 9848598,
                 gene_model = gene_anno, 
                 coaccess_cutoff = .25, 
                 connection_width = .5, 
                 collapseTranscripts = "longest" )

假如之前用的全基因组来算connections，那么这一步也会很久才能出图，所以建议先自己提取部分数据，这样子作图的速度就会快很多(因为源代码没有优化好矩阵过滤)

conns_sub <- conns[conns_matrix[,1] == "chr2" & 
                     as.numeric(conns_matrix[,2]) > 8000000 &
                     as.numeric(conns_matrix[,2]) < 10000000 
                     , ]

plot_connections(conns_sub, "chr2", 9773451, 9848598,
                 gene_model = gene_anno, 
                 coaccess_cutoff = .25, 
                 connection_width = .5, 
                 collapseTranscripts = "longest" )

connections

计算基因活跃度得分

在启动子区的开放度并不是基因表达很好的预测子。使用Cicero links，我们就能得到启动子和它相关的远距位点的总体开放情况。这种区域间开放度联合得分和基因表达的相关度更高。作者称之为基因活跃度得分（Cicero gene activity score)。可以用两个函数进行计算

build_gene_activity_matrix函数接受CDS和Cicero connection list作为输入，返回一个标准化的基因得分。注意， CDS必须在fData 表中有一列”gene”, 如果该peak在启动子区就用对应的基因名进行注释, 否则为NA. 我们可以用annotate_cds_by_site进行注释。

build_gene_activitity_matrix是未标准化的数据。需要用normalize_gene_activities进行标准化。如果你想比较不同数据集中部分数据的基因活跃度得分, 那么所有的gene activity 部分数据都需要一起标准化。如果只是想标准化其中一个数据集，那么传递一个结果即可。normalize_gene_activities也需要一个每个细胞所有开放位点的命名向量。这存放在CDS对象中的pData，在num_genes_expressed列。

可能还是不是明白，看代码吧。

第一步，是用annotate_cds_by_site对CDS对象进行注释在启动子区域的peak。

我们先从GFF文件中提取启动子的位置，也就是每个转录本的第一个外显子

# 分别提取正链和负链的注释, 仅保留第一个外显子
pos <- subset(gene_anno, strand == "+")
pos <- pos[order(pos$start),] 

# remove all but the first exons per transcript
pos <- pos[!duplicated(pos$transcript),] 
# make a 1 base pair marker of the TSS
pos$end <- pos$start + 1 

neg <- subset(gene_anno, strand == "-")
neg <- neg[order(neg$start, decreasing = TRUE),] 
# remove all but the first exons per transcript
neg <- neg[!duplicated(neg$transcript),] 
neg$start <- neg$end - 1

# 合并正链和负链
gene_annotation_sub <- rbind(pos, neg)

# 只需要4列信息
gene_annotation_sub <- gene_annotation_sub[,c("chromosome", "start", "end", "symbol")]

# 第四列命名为gene
names(gene_annotation_sub)[4] <- "gene"

用annotate_cds_by_site对cds按照gene_annotation_sub注释。会多处两列，overlap和gene。

1 2	input_cds <- annotate_cds_by_site(input_cds, gene_annotation_sub) tail(fData(input_cds))

第二步: 计算未标准化的基因活跃度得分

# generate unnormalized gene activity matrix
unnorm_ga <- build_gene_activity_matrix(input_cds, conns)

# 过滤全为0的行和列
unnorm_ga <- unnorm_ga[!Matrix::rowSums(unnorm_ga) == 0, 
                       !Matrix::colSums(unnorm_ga) == 0]

第三步: 标准化基因活跃度得分

# 你需要一个命名列表
num_genes <- pData(input_cds)$num_genes_expressed
names(num_genes) <- row.names(pData(input_cds))

# 标准化
cicero_gene_activities <- normalize_gene_activities(unnorm_ga, num_genes)

如果你有多个数据集需要进行标注化，那你就需要传递一个列表给normalize_gene_activities, 并且保证num_genes里包括所有数据集中的所有样本。

# if you had two datasets to normalize, you would pass both:
# num_genes should then include all cells from both sets
unnorm_ga2 <- unnorm_ga
cicero_gene_activities <- normalize_gene_activities(list(unnorm_ga, unnorm_ga2), 
                                                    num_genes)

标注化后的基因活跃度得分在0到1之间,

第四步: 在UMAP上可视化基因得分。由于值比较小，要先扩大1000000倍，即log2(GA * 1000000 + 1)，这样子在可视化的时候比较好看。

plot_umap_ga <- function(cds,
                         cicero_gene_activities,
                         gene,
                         dotSize = 1,
                         log = TRUE){
  #umap_df <- as.data.frame(colData(cds))
  umap_df <- as.data.frame(reducedDims(input_cds)$UMAP)
  colnames(umap_df) <- c("UMAP1","UMAP2")
  gene_pos <- which(row.names(cicero_gene_activities) %in% gene)

  umap_df$score <- log2(cicero_gene_activities[gene_pos,] * 1000000 + 1)
  if (log){
    umap_df$score <- log(umap_df$score + 1)
  }

  p <- ggplot(umap_df, aes(x=UMAP1, y=UMAP2)) +
    geom_point(aes(color=score), size = dotSize) + scale_color_viridis(option="magma") +
    theme_bw()
  return(p)
}

拟时序分析

这部分的分析其实比较简单，因为仅仅是将feature从mRNA变成了peak而已。步骤如下

数据预处理
数据降维
细胞聚类
轨迹图学习
在拟时间上对细胞进行排序

代码如下

# 读取数据
temp <- textConnection(readLines(gzcon(url("http://staff.washington.edu/hpliner/data/kidney_data.txt.gz"))))set.seed(2017)
input_cds <- estimate_size_factors(input_cds)
#1 数据预处理
input_cds <- preprocess_cds(input_cds, method = "LSI")
#2 降维
input_cds <- reduce_dimension(input_cds, reduction_method = 'UMAP', 
preprocess_method = "LSI")
#3 细胞聚类
input_cds <- cluster_cells(input_cds)
#4 轨迹图学习
input_cds <- learn_graph(input_cds)
#5 对细胞进行排序
input_cds <- order_cells(input_cds, 
                         root_cells = "GAGATTCCAGTTGAATCACTCCATCGAGATAGAGGC")

最后可视化

1	plot_cells(input_cds, color_cells_by = "pseudotime")

差异开放分析

当你的细胞在拟时间中排列之后，你可能会想知道基因组哪些区域会随着时间发生变化。

如果你已经知道哪些区域会随着时间发生变化，你可以用plot_accessibility_in_pseudotime进行可视化

input_cds_lin <- input_cds[,is.finite(pseudotime(input_cds))]
plot_accessibility_in_pseudotime(input_cds_lin[c("chr1_3238849_3239700", 
                                                 "chr1_3406155_3407044", 
                                                 "chr1_3397204_3397842")])

当然我们也可以从头分析. 先用aggregate_by_cell_bin对细胞根据相似度进行聚合

# First, assign a column in the pData table to umap pseudotime
pData(input_cds_lin)$Pseudotime <- pseudotime(input_cds_lin)
pData(input_cds_lin)$cell_subtype <- cut(pseudotime(input_cds_lin), 10)
binned_input_lin <- aggregate_by_cell_bin(input_cds_lin, "cell_subtype")

之后用Monocle3的fit_models函数寻找差异开发的区域。

# For speed, run fit_models on 1000 randomly chosen genes
set.seed(1000)
acc_fits <- fit_models(binned_input_lin[sample(1:nrow(fData(binned_input_lin)), 1000),], 
                       model_formula_str = "~Pseudotime + num_genes_expressed" )
fit_coefs <- coefficient_table(acc_fits)

# Subset out the differentially accessible sites with respect to Pseudotime
pseudotime_terms <- subset(fit_coefs, term == "Pseudotime" & q_value < .05)
head(pseudotime_terms)

参考资料

官方教程:https://cole-trapnell-lab.github.io/cicero-release/docs_m3

2019-08-17

生信软件工具箱

motifmatchr:在R语言中分析peak中里是否有motif匹配

motifmatchr的作用就是分析众多的序列和众多的motifs, 从中找到哪个序列包含哪个motif. 它的核心函数就是matchMotifs，最大特点就是快，因为它用的是MOODS C++库用于motif匹配。

尽管Bioconductor上也有很多工具能够做motif匹配，比如说Biostrings::mathcPWM, TFBSTools::searchSeq，但是motifmatchr更适合于分析许多不同的序列包含许多不同的motif。例如，当分析ChIP-seq或者ATAC-seq数据时, 你可能会想知道在哪个peak里有哪种类型的motif.

R包安装和加载

1 2	library(motifmatchr) library(GenomicRanges)

matchMotifs

motifmatchr的核心函数是matchMotifs，因此了解这个函数的输入数据是什么格式就行了。必须的两个输入是

位置权重矩阵(position weight matrices, PWM)或位置频率矩阵(position frequency matrices, PFM), 保存在PWMatrix, PFMatrix, PWMatrixList或PFMatrixList
一组基因组范围(GenomicRanges或RangedSUmmarizedExperiment对象)或一组序列(DNAStringSet, DNAString 或简单的字符串向量)

PWM或PFM矩阵

Motif的PWM或PFM矩阵可以从JASPAR, CIS-BP下载。

例如，我从CIS-BP下载拟南芥Arabidopsis_thaliana_2019_08_17_2_32_am.zip，解压缩后里有一个pwm_all_motifs文件夹，里面的文本文件存放的就是我所需要的PWM矩阵, 下一步根据这些文件构建出matchMotifs的输入对象

motif_dir <- "Arabidopsis_thaliana_2019_08_17_2_32_am/pwms_all_motifs"
PWList <- list()
for (file in list.files(motif_dir, pattern = ".txt")){
  df <- read.table(file.path(motif_dir, file), 
                   header = T,
                   row.names = 1)
  mt <- as.matrix(df)
  if (nrow(mt) ==0) next
  motif_id <- substr(file, 1,6)
  PWList[[motif_id]] <- PWMatrix(ID = motif_id, profileMatrix = t(mt))
}

PWMatrixLists <- do.call(PWMatrixList,PWList)

对于JASPAR，我们有更加方便的提取方法

library(JASPAR2018)
species <- "Arabidopsis thaliana"
collection <- "CORE"
opts <- list()
opts["species"] <- species
opts["collection"] <- collection
out  <- TFBSTools::getMatrixSet(JASPAR2018, opts)

if (!isTRUE(all.equal(TFBSTools::name(out), names(out)))) 
  names(out) <- paste(names(out), TFBSTools::name(out), 
                      sep = "_")
motif <- out

JASPAR的motif比较可靠，因此motif相对比较少。

给定范围或序列

这部分的输入就比较容易获取，你可以提供MACS2的输出BED，利用rtracklayer::import.bed()读取构建成一个GRanges对象。因为是提供的GRanges对象，那么还需要额外设置一个参数genome, 利用Biostrings::readDNAStringSet()读取一个参考基因组序列就行了。

或者用bedtools先根据bed文件提取数据，然后利用Biostrings::readDNAStringSet()读取

示例数据

我们以包中提供的数据为例，进行演示

加载实例的motif和GRanges对象

# load some example motifs
data(example_motifs, package = "motifmatchr") 

# Make a set of peaks
peaks <- GRanges(seqnames = c("chr1","chr2","chr2"),
                 ranges = IRanges(start = c(76585873,42772928,100183786),
                                  width = 500))

获取motif在peak中的位置

1	motif_ix <- matchMotifs(example_motifs, peaks, genome = "hg19")

用motifMatches函数提取匹配矩阵

1	motifMatches(motif_ix)

输出结果是一个稀疏矩阵

3 x 3 sparse Matrix of class "lgCMatrix"
     MA0599.1_KLF5 MA0107.1_RELA MA0137.3_STAT1
[1,]             |             |              .
[2,]             |             .              |
[3,]             |             .              |

其中的.就表示存在motif匹配。

我们还可以提取motif在peak中的位置

1
2
3

# Get motif positions within peaks for example motifs in peaks 
motif_ix <- matchMotifs(example_motifs, peaks, genome = "hg19",
                         out = "positions")

其他参数

除了必须的motif信息和你的序列信息输入外，还有一些其他的参数可以做一些设置。

背景核苷酸频率, 默认是bg=subject, 也就是你的输入数据作为背景，也可设置成genome或even
P值: 用于判断匹配是否足够可信的参数，默认是0.00005，基本上不用修改
输出信息: matchMotifs可以返回三种可能输出，matches, scores 和 positions

总的来说，这个R包是一个比较简单的工具，比较容易上手使用。

2019-08-17

R语言的稀疏矩阵学习记录

一个很大的矩阵， 320127 行, 8189列，假如用一个全为0的普通矩阵来存储，需要用到9.8Gb

cols <- 8189
rows <- 320127
mat <- matrix(data = 0, nrow=320127, ncol = 8189)
print(object.size(mat), unit="GB")
# 19.5 Gb
mat <- matrix(data = 0L, nrow=320127, ncol = 8189)
print(object.size(mat), unit="GB")
# 9.8 Gb这里的0其实也要区分

这里的0L表示数据类型是integer，默认是numeric. 这两者最大的区别在于，当你用320127L * 8189L，你会得到一个NA，而320127 * 8189不会

如果用稀疏矩阵保存的话

mat <- Matrix(data = 0L, nrow=320127, ncol = 8189, sparse = TRUE)
print(object.size(mat), unit="GB")
#0 Gb
dim(mat)
#[1] 320127   8189

虽然行列数一样，但是稀疏矩阵几乎不占用任何内存。而且普通矩阵支持的运算，比如说求行和，求列和，提取元素的操作，在稀疏矩阵矩阵也是可以的，只不过会多花一点点时间而已。同时还有很多R包支持稀疏矩阵，比如说glmnet，一个做lasso回归的R包。

虽然看起来稀疏矩阵很美好，但是在R语言中那么大的稀疏矩阵的部分操作会出错

1
2
3

> mat2 <- mat + 1
Error in asMethod(object) : 
  Cholmod error 'problem too large' at file ../Core/cholmod_dense.c, line 105

即便是我想把它用as.matrix转回普通矩阵，它也报错了

1
2
3

> mat3 <- Matrix::as.matrix(mat)
Error in asMethod(object) : 
  Cholmod error 'problem too large' at file ../Core/cholmod_dense.c, line 105

既然现成的as.matrix无法处理，那怎么办呢？最简单粗暴的方法就是新建一个普通矩阵，然后对稀疏矩阵进行遍历，将稀疏矩阵的值挨个放回到的普通矩阵上。

mat2 <- matrix(data = 0, nrow=320127, ncol = 8189)
for (i in seq_len(nrow(mat))){
    for (j in seq_len(ncol(mat))){
        mat2[i,j] <- mat[i,j]
    }
}

那么这大概要多少时间呢？反正我的电脑跑了2个小时也没有跑完，所以你也别测试了。

那有没有办法可以加速呢？加速的方法就是减少for循环的次数，因为我们是一个稀疏矩阵，大部分的空间都是0，我们只需要将不为0的部分赋值给新矩阵即可。

这需要我们去了解下稀疏矩阵的数据结构

> str(mat)
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  ..@ i       : int(0) 
  ..@ p       : int [1:8190] 0 0 0 0 0 0 0 0 0 0 ...
  ..@ Dim     : int [1:2] 320127 8189
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
  ..@ x       : num(0) 
  ..@ factors : list()

@Dim记录矩阵的维度信息, @Dimnames记录行名和列名, @x记录不为0的数值。@i记录不为0的行索引，和@x对应，这里全为0，所以不记录。@p比较复杂，并不是简单的记录不为0值的列索引，看文档也不知道是啥，不过通过检索可以找到它和不为0值的列索引的换算关系。

因此代码优化为

row_pos <- mat@i+1
col_pos <- findInterval(seq(mat@x)-1,mat@p[-1])+1
val <- mat@x
    
for (i in seq_along(val)){
    tmp[row_pos[i],col_pos[i]] <- val[i]
}

可以将其封装为一个函数

as_matrix <- function(mat){

  tmp <- matrix(data=0L, nrow = mat@Dim[1], ncol = mat@Dim[2])
  
  row_pos <- mat@i+1
  col_pos <- findInterval(seq(mat@x)-1,mat@p[-1])+1
  val <- mat@x
    
  for (i in seq_along(val)){
      tmp[row_pos[i],col_pos[i]] <- val[i]
  }
    
  row.names(tmp) <- mat@Dimnames[[1]]
  colnames(tmp) <- mat@Dimnames[[2]]
  return(tmp)
}

如果速度还需要提高，那么可能就需要Rcpp上场了. 我参考着http://adv-r.had.co.nz/Rcpp.html写了一个简单的代码

Rcpp::sourceCpp(code='
#include <Rcpp.h>
using namespace Rcpp;


// [[Rcpp::export]]
IntegerMatrix asMatrix(NumericVector rp,
                       NumericVector cp,
                       NumericVector z,
                       int nrows,
                       int ncols){

  int k = z.size() ;

  IntegerMatrix  mat(nrows, ncols);

  for (int i = 0; i < k; i++){
      mat(rp[i],cp[i]) = z[i];
  }

  return mat;
}
' )


as_matrix <- function(mat){

  row_pos <- mat@i
  col_pos <- findInterval(seq(mat@x)-1,mat@p[-1])
  
  tmp <- asMatrix(rp = row_pos, cp = col_pos, z = mat@x,
                 nrows =  mat@Dim[1], ncols = mat@Dim[2])

  row.names(tmp) <- mat@Dimnames[[1]]
  colnames(tmp) <- mat@Dimnames[[2]]
  return(tmp)
}

如果之前的矩阵有78945836个元素，system.time显示只需要40s。

2019-08-16

生信软件工具箱

chromVAR:从单细胞表观数据中推断转录因子相关开放区域

chromVAR是一个用于分析稀疏染色质开放的R包。chromVAR的输入文件包括，ATAC-seq处理后的fragments文件（过滤重复和低质量数据）, DNAse-seq实验结果，以及基因组注释（例如motif位置）

chromVAR先根据所有细胞或者样本的平均情况来计算期望开放性，然后用它来计算每个注释，每个细胞或样本的偏差，最后对开放进行纠正。

安装加载

安装R包

BiocManager::install("GreenleafLab/chromVAR")
BiocManager::install("motifmatchr")
BiocManager::install(SummarizedExperiment)
BiocManager::install(BiocParallel)
BiocManager::install("JASPAR2018") # JASPAR 2018数据库

chromVAR的运行需要先加载下面这些R包

library(chromVAR)
library(motifmatchr)
library(Matrix)
library(SummarizedExperiment)
library(BiocParallel)
set.seed(2019)

chromVAR能够使用多核进行并行运算，调用方法如下

# 全部用户
register(MulticoreParam(8, progressbar = TRUE))
# Windows用户，MulticoreParam不好用就用SnowParam
register(SnowParam(workers = 1, type = "SOCK"))

注意: 不运行的话，后续代码可能会报错

读取输入

chromVAR接受的输入是落入开放区域的read数统计表。有许多软件可以做到，chromVAR也提供了相应的方法。

首先要提供一个peak文件，文档建议这个peak文件存放的peak为等宽非重叠，建议宽度在250-500 bp之间。peak文件可以利用已有的bulk ATAC-seq或DNAse-seq数据来获取。对于来自于多个样本的peak，需要先用filterPeaks函数保证peak之间不重合。

使用getPeaks读取peak

1 2	peakfile <- system.file("extdata/test_bed.txt", package = "chromVAR") peaks <- getPeaks(peakfile, sort_peaks = TRUE)

MACS2分析结果里会提供narrowpeak格式文件，chromVAR提供了readNarrowpeaks函数进行读取。

随后用getCounts函数基于BAM文件和加载的peak获取count

bamfile <- system.file("extdata/test_RG.bam", package = "chromVAR")
fragment_counts <- getCounts(bamfile, 
                             peaks, 
                             paired =  TRUE, 
                             by_rg = TRUE, 
                             format = "bam", 
                             colData = DataFrame(celltype = "GM"))

bamfile可以有多个bam文件路径，bam文件的类型和colData对应。此外by_rg定义是否要根据BAM文件中的RG对输入分组。

实际演示的时候，我们官方提供的示例数据

1	data(example_counts, package = "chromVAR")

因为这是一个*SummarizedExperiment 对象，因此我们就能用assay,colData和rowData 了解这个数据集

主体是一个存放count的dgCMatrix

assay(example_counts)[1:5,1:2]

5 x 2 sparse Matrix of class "dgCMatrix"
     singles-GM12878-140905-1 singles-GM12878-140905-2
[1,]                        .                        1
[2,]                        .                        .
[3,]                        .                        .
[4,]                        .                        .
[5,]                        .                        .

行是特征(feature), 这里就是peak

head(rowData(example_counts), n=2)

DataFrame with 2 rows and 3 columns
      score      qval        name
  <integer> <numeric> <character>
1       259  25.99629   GM_peak_6
2        82   8.21494   H1_peak_

列表示样本，

head(colData(example_counts), n=2)

DataFrame with 2 rows and 2 columns
                           Cell_Type     depth
                         <character> <integer>
singles-GM12878-140905-1          GM      9220
singles-GM12878-140905-2          GM      9401

前期准备

分析peak中的GC含量

GC含量信息用于确定哪些peak可能是背景， addGCBias返回更新后的SummarizedExperiment 。其中genome支持BSgenome, FaFile和DNAStringSet对象。

library(BSgenome.Hsapiens.UCSC.hg19)
example_counts <- addGCBias(example_counts, 
                            genome = BSgenome.Hsapiens.UCSC.hg19)
head(rowData(example_counts))

过滤输入(单细胞)

如果是处理单细胞数据，建议过滤测序深度不够，或者是信噪比低不够(也就是peak中所占read比例低)的数据，主要靠两个参数,min_in_peaks和min_depth.

1 2	counts_filtered <- filterSamples(example_counts, min_depth = 1500, min_in_peaks = 0.15, shiny = FALSE)

然后我们用filterSamplesPlot对过滤情况进行可视化

1
2
3

filtering_plot <- filterSamplesPlot(example_counts, min_depth = 1500, 
                                    min_in_peaks = 0.15, use_plotly = FALSE)
filtering_plot

确认标准后，就可以过滤

1	counts_filtered <- filterPeaks(counts_filtered, non_overlapping = TRUE)

获取Motifs和分析包含motif的peak

可以用getJasparMotifs在JASPAR数据库中提取motif信息，但是其中getJasparMotifs只是TFBSTools::getMatrixSet一个简单的封装而已, 默认用的是JASPAR2016, 建议通过下面的代码获取最新版本的数据。

# BiocManager::install("JASPAR2018")
library(JASPAR2018)
species <- "Homo sapiens"
collection <- "CORE"
opts <- list()
opts["species"] <- species
opts["collection"] <- collection
out  <- TFBSTools::getMatrixSet(JASPAR2018, opts)

if (!isTRUE(all.equal(TFBSTools::name(out), names(out)))) 
  names(out) <- paste(names(out), TFBSTools::name(out), 
                      sep = "_")
motif <- out

collection参数中接受: “CORE”, “CNE”, “PHYLOFACTS”, “SPLICE”, “POLII”, “FAM”, “PBM”, “PBM_HOMEO”, “PBM_HLH” 选项。

获取Motif的另外一种方式是用chromVARmotifs, 这个R包的主要功能就是整合了一些motif, 通过devtools::install_github("GreenleafLab/chromVARmotifs")进行安装。

之后用motifmatchr中的macthMotifs去分析peak中motif包含情况。默认会返回一个SummarizedExperiment 对象，就是一个稀疏矩阵，来提示是不是匹配了motif

1
2
3

library(motifmatchr)
motif_ix <- matchMotifs(motifs, counts_filtered, 
                        genome = BSgenome.Hsapiens.UCSC.hg19)

关键参数是p.cutoff用于设置严格度，默认是0.00005其实能够返回比较合理的结果。如果需要一些额外信息，可以通过参数out来调整，例如out=positions表示返回实际的匹配位置

偏离度和变异度分析(Deviations and Variability )

上面都是准备阶段，计算deviations和Variability才是这个软件的重点。

计算偏离度

函数computeDeviations返回的SummarizedExperiment 对象, 包含两个”assays”,

偏差纠正后的开放偏离度:assays(dev)$deviations
偏离度的Z-score: assays(deviations)$z: 考虑到了当从基因组上随机抽取一些GC含量类似的片段分析时，有多大概率会得到该结果。

1	dev <- computeDeviations(object = counts_filtered, annotations = motif_ix)

compuateDeviations支持背景peak的输入，用于的偏移度得分进行标准化，默认会自动计算不会返回结果。你可以选择用getBackgroundPeaks分析背景，然后传递给computeDeviations

1
2
3

bg <- getBackgroundPeaks(object = counts_filtered)
dev <- computeDeviations(object = counts_filtered, annotations = motif_ix,
                         background_peaks = bg)

变异度

computeVariability会返回一个数据框data.frame，里面有变异度, 该变异度的bootstrap置信区间和衡量拒绝原假设的概率(即 > 1)

1 2	variability <- computeVariability(dev) plotVariability(variability, use_plotly = FALSE)

偏离度可视化

可以利用tSNE可视化偏离度

tsne_results <- deviationsTsne(dev, threshold = 1.5, perplexity = 10)
tsne_plots <- plotDeviationsTsne(dev, tsne_results, 
                                 annotation_name = "TEAD3", 
                                 sample_column = "Cell_Type", 
                                 shiny = FALSE)
cowplot::plot_grid(tsne_plots[[1]],tsne_plots[[1]] )

参考资料

数据库JASPAR是一个开放性数据库，用于存放人工审查的非冗余转录因子(TF)结合谱，数据格式为PFM(position frequency matrices ), TFFM( TF flexible models)
chromVAR官方文档
chromVAR: inferring transcription-factor-associated accessibility from single-cell epigenomic data

2019-08-15

生物信息学格式和转换工具整理

生信会用到很多格式，不同软件会要求不同的输入格式，所以得要自己整理下。

FASTA和FASTQ

SAM和BAM

GTF和GFF

BED

BedGraph

Wig和BigWig

WIG（wiggle）设计的目的是在基因组浏览器上展示连续性数据。Wiggle数据要求元素等长。对于稀疏或者元素之间不等长，建议用BedGraph格式。而BigWig则是二进制版，体积小，适合传输。

举个例子

variableStep chrom=chr2
300701  12.5
300702  12.5
300703  12.5
300704  12.5
300705  12.5

和

1 2	variableStep chrom=chr2 span=5 300701 12.5

在基因组浏览器上展示时都是2号染色体的300701-300705的值为12.5