xuzhougeng blog

因为有些时候上传的是BAM文件

在做基因组注释时，有些时候需要将近源物种的蛋白序列回帖到基因组上，之前用到GeneWise或者Exonerate，但是这两个工具的岁数都太大了，且速度也不够快。好在，李恒团队出手了，开发了miniprot，从命名上，就知道他会跟minimap一样高效，好用。这里的ref.fna是参考基因组序列。之所以要这样子做，是因为miniprot需要将基因组翻译成氨基酸序列，这一部分比较费时间。输出的gff文件就可以直接用IGV进行可视化，或者用其他软件解析，作为下游分析的输入。安装非常简单，仅依赖于zlib。

通过后续的分析发现，只需要保证输入的fasta是多行形式，而不是单行形式，那么原来的流程的分析 速度并不会受到影响。可以使用 seqkit seq -w 80 input.fasta > output.fasta 的方式来调整输入fasta的格式。

2019年发表在Nature上的文章【The single-cell transcriptional landscape of mammalian organogenesis】在方法部分提到，使用NNLS(non-negative linear-square)回归的方法分析两个细胞图谱数据集中相关细胞类型。这个方法，在我搜索的中文教程中都没有出现过，所以这里以两个pbmc的数据集为例进行介绍，如何复现文章的方法。10x的细胞数据集的预处理部分不做过多介绍， 如下代码以10x官网提供的数据为例libra

代码和数据都在GitHub上，见 https://github.com/xuzhougeng/rice_annotaiton

在如何在IGV上使用BLAT搜索非模式物种中，我们讨论如何用Apache的CGI服务器，来响应IGV发送的BLAT请求。考虑到我们大部分的时候都是个人使用，并不需要Apache这种重量级的Web服务器，因此完全可以省去这一组件。我之前学过Python和Django和Flask这两种网页应用框架，其中Flask比较轻量，非常适合我们这种小型应用。因此，我就用Flask编写了一个Blat网页应用，用来响应IGV的请求。首先，安装Flask（为了避免冲突，我用了虚拟环境)mkdir myprojectcd

IGV提供了BLAT，用于进行序列搜索，但可惜我一直用不上，因为它默认是调用了UCSC的CGI工具，将我们的输入序列发送到https://genome.ucsc.edu/cgi-bin/hgBlat处理，处理后返回JSON文件用于展示。因此，除非我们自己搭建一个UCSC类似的网站，否则，无法用到IGV的这个功能。我觉得肯定不只是我一个人有这个问题，所以我就去谷歌上用关键词 “custom genome BLAT IGV” 进行检索, 果然发现有很多人都有类似的需求，我找到最早一条是2016年，但是距今.

hifiasm是一个能有效利用PacBio HiFi测序技术，在分型组装图(pahsed assembly gprah)中可靠的表示单倍体信息的算法。流程介绍hifiasm的分析流程如下，主要分为3个阶段第一阶段：通过所有序列的相互比对，对前在测序错误进行纠正。如果一个位置只存在两种碱基类型，且每个碱基类型至少有3条read支持，那么这个位置会被当作杂合位点，否则，视作测序错误，将被纠正。第二阶段：根据序列之间的重叠关系，构建分型的字符串图(phased string graph)。其中调整朝向的序

前些天发现自己系统里的一些R包实在太久了，于是头脑一热就把所有包都升级了。结果运行Seurat的时候，就卡在 FindNeighbor 了。运行到一半就提示如下报错Error in validObject(.Object) : invalid class “Graph” object: superclass "Mnumeric" not defined in the environment of the object's class于是我发动科研工作者的技能，search，查了查相关的资料，发

这部分内容是Extending ggplot2的学习笔记，大部分内容都是原文的简单翻译。所有的ggplot2对象都建立自"ggproto"这套面向对象编程系统，因此想要创建出自己的一套图层，而不是简单的对已有图层进行累加，那么就需要学习"ggproto"。创建新的stat最简单的stat我们会从一个最简单的stat开始: 根据已有的一组点，用一个凸壳(convex hull)包围他。第...

外周血单核细胞(Peripheral blood monoculear cell, PBMC), 包括淋巴细胞(T细胞，B细胞和自然杀伤(NK)细胞)和单核细胞。而红细胞和血小板没有细胞核，而粒细胞(granulocytes) 包括中性粒细胞、嗜碱性粒细胞和嗜酸性粒细胞，有多叶核，所以不包括在PBMC中。人类的PBMC主要由淋巴细胞组成，之后是单核细胞以及一小部分的树突状细胞(dendritic...

inferCNV用与探索肿瘤单细胞RNA-seq数据，分析其中的体细胞大规模染色体拷贝数变化(copy number alterations, CNA), 例如整条染色体或大片段染色体的增加或丢失(gain or deletions)。工作原理是，以一组"正常"细胞作为参考，分析肿瘤基因组上各个位置的基因表达量强度变化. 通过热图的形式展示每条染色体上的基因相对表达量，相对于正常细胞，肿瘤基因组...

我的生物信息学技能基本上都是自学的，在刚开始的阶段我也和大部分的人一样，会遇到很多问题，然后想着去一些群里提问希望得到解答，当然结果也是和大部分人的一样，问题不会得到答案。后来，我通过看书和搜索慢慢的有了一些基础，通过项目实践对生信有了一点了解，加上自己也喜欢把自己的学习笔记公开出来，就逐渐有人开始向我提问。有些问题我会想去回答，而有些问题会被我直接忽视。如同当年我提出的问题一样，有些能够有人...

mclust(Model-based clustering) 能够基于高斯有限混合模型进行聚类，分类以及密度估计(density estimation)。对于具有各种协方差结构的高斯混合模型，它提供了根据EM算法的参数预测函数。它也提供了根据模型进行模拟的函数。还提供了一类函数，整合了基于模型的层次聚类，混合估计的EM算法，用于聚类、密度估计和判别分析中综合性策略的贝叶斯信息判别标准。最后还有一...

在学习「数据挖掘导论」的数据预处理时，里面谈到了变量变换，我联想到了在基因表达量分析时的常见操作，例如FPKM，TPM，CPM，log对数变换。 比如说在文章里面会见到如下的描述The size factor of each cell was computed using a pooling strategy implemented in the R function computeSumFa...

使用ALLHiC基于HiC数据辅助基因组组装基因组组装大致可以分为三步（1）根据序列之间的重叠情况构建出contig，（2）基于二代的mate pair文库或光学图谱将contig搭建成scaffold，（3）对scaffold进行排序和调整方向得到最终的准染色体级别的基因组。目前的三代测序组装能够搞定第一步和第二步。而在将contig/scaffold提升至准染色体水平上，有4种方案可选。一...

ALLHiC续: 如何构建Allele.ctg.tableAllele.ctg.table是ALLHiC用于过滤多倍体基因组中因等位序列相似引起的HiC噪音的必要输入。构建的方法有很多种，这里列举两个方法。方法一: 基于BLAST基于BLAST的方法需要你提供一个染色体级别组装的近缘物种。需要提供4个文件，近缘物种的CDS和GFF文件，多倍体物种的CDS和GFF文件。最重要的是，你得保证C...

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

实用Seurat自带的热图函数DoHeatmap绘制的热图，感觉有点不上档次，于是我尝试使用ComplexHeatmap这个R包来对结果进行展示。个人觉得好的热图有三个要素聚类: 能够让别人一眼就看到模式注释: 附加注释能提供更多信息配色: 要符合直觉，比如说大部分都会认为红色是高表达，蓝色是低表达在正式开始之前，我们需要先获取一下pbmc的数据，Seurat提供了R包Seura...

Tribonacci序列Tn定义:T0=0, T1=1, T2=1, n>=0时，Tn 3 = Tn Tn 1 Tn 2限制条件是: 0

今天，我开发的R包，binmapr，终于上线CRAN了，也就意味着你可以通过install.packages("binmarp")的方式进行安装了。不过这个R包目前估计也就只有我会用，所以近期会根据几个具体案例来介绍如何使用这个R包。这篇讲讲我写这个R包的一些心得体会写函数要尽早这句话是听Y叔说的。当你在写代码的时候，发现一个操作要反复进行，那么最好的方式就是把它封装成一个函数，放在一个scr...

SCALE全称是Single-Cell ATAC-seq analysis vie Latent feature Extraction, 从名字中就能知道这个软件是通过隐特征提取的方式分析单细胞ATAC-seq数据。在文章中，作者从开发者的角度列出了目前的scATAC-seq分析软件，chromVAR, scABC, cisTopic, scVI，发现每个软件都有一定的不足之处，而从我们软件使用...

「文献」杂合基因组的策略思路文献出处: Sequencing a Juglans regia × J. microcarpa hybrid yields high-quality genome assemblies of parental specie文章的亮点在于通过对一个F1子代进行三代测序，之后利用BioNano组装出两个亲本的光学图谱，最后根据光学图谱从F1中拆分出两套单倍型。杂合基...

服务器上安装RStudio-server如果想在服务器上安装一个RStudio-server，你需要先保证自己拥有管理员权限，之后参考如何在服务器上安装最新的R安装R语言，一定要注意在./configure的时候加上--enable-R-shlib，否则后续会出错。RStudio-server分为两种版本，一种是开源免费版，另一个是商业专业版本。个人觉得两者最大的区别在于，商业版支持在多个版本...

## 目前已有的框架A review of bioinformatics pipeline framework 的作者对已有的工具进行很好的分类作者的看法：implicit，也就是Make rule语法更适合用于整合不同执行工具基于配置的流程更加稳定，也比较适合用于集群分配任务。最后作者建议是：如果实验室既不是纯粹的生物学试验（不需要workbench这种UI界面），也不需要高...

尽管目前已测序的物种已经很多了，但是对于一些动辄几个G的复杂基因组，还没有某个课题组有那么大的经费去测序，所以仍旧缺少高质量的完整基因组，那么这个时候一个高质量的转录组还是能够凑合用的。二代测序的组装结果只能是差强人意，最好的结果就是不要组装，直把原模原样的转录组给你是最好的。PacBio Iso-Seq 做的就是这件事情。只不过Iso-Seq测得是转录本，由于有些基因存在可变剪切现象，所以所有...

这里的新版指的是PacBio公司在2018年9月发布pb-assembly, 而这篇文章是在2018年9月30日发的。今年早些时候在参加三代培训时，听说PacBio会在今年对Falcon进行一些改变。前几天我在读 readthedocs上的Falcon文档时，发现了文档页面上方出现了这样两栏提醒其中pb_assembly就是新的FALCON组装套装在GitHub上的使用文档，经过了几天的...

基因组survey在组装基因组之前一定要先对要组装的物种有一个大致的了解，判断其复杂程度, 标准如下简单基因组: 杂合度低于0.5%, GC含量在35%~65%, 重复序列低于50%二倍体普通基因组: 杂合度在0.5%~1.2%中间，重复序列低于50%。或杂合度低于0.5%，重复序列低于65%高复杂基因组: 杂合度>1.2% 或 重复率大于65%k-mers最简单...

基因组组装基因组组装一般分为三个层次，contig, scaffold和chromosomes. contig表示从大规模测序得到的短读(reads)中找到的一致性序列。组装的第一步就是从短片段(pair-end)文库中组装出contig。进一步基于不同长度的大片段(mate-pair)文库，将原本孤立的contig按序前后连接，其中会调整contig方向以及contig可能会存在开口(gap...

我们用来练手的文章发表在 Nature Communication ，”High contiguity Arabidopsis thaliana genome assembly with a single nanopore flow cell”, 非常不要脸的说，这篇文章是我师爷实验室发的。简单讲讲故事内容，就是他们实验室买了一台nanopore仪器，就是下面这台， 目前仪器价格国内是8K左右...

Basic local alignment search tool (BLAST)包括：blastn, blastp, blastx, tblastn, tblastx等. 使用conda安装即可。conda install -c bioconda blast# blast安装perl模块的方法conda isntall perl-digest-md5BLAST的主要理念...

基因组组装完成后，或者是完成了草图，就不可避免遇到一个问题，需要对基因组序列进行注释。从头注释(de novo prediction)：通过已有的概率模型来预测基因结构，在预测剪切位点和UTR区准确性较低同源预测(homology-based prediction)：有一些基因蛋白在相近物种间的保守型搞，所以可以使用已有的高质量近缘物种注释信息通过序...

如何使用fastq-dump转换SRA格式做生信的基本上都跟NCBI-SRA打过交道,尤其是fastq-dump大家肯定不陌生.NCBI的fastq-dump软件一直被大家归为目前网上文档做的最差的软件之一”,而我用默认参数到现在基本也没有出现过什么问题,感觉好像也没有啥问题, 直到今天看到如下内容, 并且用谷歌搜索的时候,才觉得大家对fastq-dump的评价非常很到位.我们一般使用...

训练 ab initio 基因预测工具(以SNAP为例)对于一个新的物种而言，你大概率是没有一个高质量的基因模型去进行基因预测。但是我们可以利用EST序列（少部分物种估计有）、二代测序数据、同源物种蛋白序列，先直接用Maker做基因注释，尽管得到的模型可能不是特别的完美，但可以作为输入反复迭代运行Maker，从而提高最终的表现。这次使用的是下载的练习数据集（见附录）cd ~...

在基因组注释上，MAKER算是一个很强大的分析流程。能够识别重复序列，将EST和蛋白序列比对到基因组，进行从头预测，并在最后整合这三个结果保证结果的可靠性。此外，MAKER还可以不断训练，最初的输出结果可以继续用作输入训练基因预测的算法，从而获取更高质量的基因模型。Maker的使用比较简单，在软件安装成后，会有一个”data”文件夹存放测试数据ls ~/opt/biosoft/ma...

最近我拿到了一批Bionano数据，用关键字 “Bionano+组装” 进行检索时，并没有发现任何的教程，所以这应是中文网络世界里第一个Bionano数据分析系列Bionano技术简单来说，就是给分子加上荧光标记，然后拍照，所以最原始的下机数据就是TIFF格式，但是用户拿到的一般都是经AutoDetect/IrysView 转换过的BNX格式。这篇文章主要就是讲讲BNX格式的具体含义。根据...

从这部分开始，就开始涉及一些软件的操作和数据分析，因此在进入正文之前，我们需要准备好环境。环境准备第一步：从 https://bionanogenomics.com/library/datasets/下载人类测试数据集，以及对应的NA12878人类基因组。wget http://bnxinstall.com/publicdatasets/DLS/20180413_NA12878_S3_co...

在2年之前我写过一篇教程介绍MUMmer软件的使用方法，可以通过如何使用MUMmer比对大片段序列阅读。MUMmer作为一个比对工具，它的主要功能就是寻找两个序列的相似之处，至于如何展示结果，并不是它的主要目标。这篇文章将会介绍如何基于MUMMER的输出结果进行可视化。首先是下载数据，我们用细菌的基因组作为案例wget http://mummer.sourceforge.net/examp...

使用BRAKER2进行基因组注释BRAKER2是一个基因组注释流程，能够组合GeneMark，AUGUSTUS和转录组数据。在使用软件之前，有几点需要注意下尽量提供高质量的基因组。目前随着三代测序价格下降，这一点问题不大。基因组命名应该简单，最好就是">contig1"或">tig000001"基因组需要屏蔽重复序列默认参数通常表现效果就很好，但是也要根据物种来一定要对注释结果进行...

「文献」多倍体植物基因组测序组装当前策略文献地址: Current Strategies of Polyploid Plant Genome Sequence Assembly基因组多倍化主要发生在被子植物中。很多多倍体植物都在农业生产上有重大的价值，例如小麦(Triticum aestivum)，花生(Arachis hypogaea)，十字花科，马铃薯(Solanum tuberosum)...