基于 genmap 计算基因组 mappability
基因组的 mappability 是指测序产生的短读长能够被唯一且准确比对到参考基因组特定区域的能力。它反映了基因组不同区域因序列特征（如重复序列、低复杂度区域等）对比对结果可靠性的影响。例如，高度重复区域（如 Alu 元件）的 mappability 通常较低，因为 reads 可能匹配到多个位置。
在实际分析中，有时我们会过滤低 mappability 区域的突变，避免因比对错误导致的假阳性误判。本文将介绍如何使用生物信息学工具 GenMap 获得基因组的 mappability 信息。
GenMapGenMap github page：
https://github.com/cpockrandt/genmap
GenMap paper：
GenMap: ultra-fast computation of genome mappability
https://doi.org/10.1093/bioinformatics/btaa222
相较于先前的 mappability 计算工具（例如 GEM），GenMap 的搜索策略有所优化并且算法不依赖于启发式方法，因此具有更高的运行效率和 ...
VCF 文件转坐标的方法汇总及注意事项
历史数据的整合分析往往涉及不同时期使用不同参考基因组版本产生的数据。例如，某些较早的研究可能基于 hg19/GRCh37 进行分析，而新的研究多采用 hg38/GRCh38。为了进行比较分析，需要将这些数据统一到同一个版本。
坐标转换是通过 chain file 来实现的，该文件记录了两个基因组版本之间的对应关系，包含了基因组重排、插入和删除等信息。转换工具会根据这些信息将原始坐标映射到目标版本的相应位置。
正文本文将介绍如何使用两款主流工具 CrossMap 和 GATK LiftoverVcf 以及最近的一个新工具 BCFtools liftover 进行 VCF 文件的坐标转换。在说明之前，首先简要介绍三个软件的差别：
CrossMap 除了 VCF 的坐标映射外，也支持很多其他文件格式的坐标转换。对于 VCF 文件而言，CrossMap 并不要求其具有完备的 Meta-information 和 Header line 部分，只需要突变数据行符合格式即可。
GATK LiftoverVcf 则主要支持对 VCF 文件的坐标进行转换。其要求 VCF 文件具有完整的 Meta-i ...
bedtools merge 不合并相邻区间（仅合并重叠区间）
关于 bedtools merge 的详细参数设定可见：
bedtools merge documentation：
https://bedtools.readthedocs.io/en/latest/content/tools/merge.html
默认情况下，bedtools merge 会对所有重叠区间及相邻区间进行合并，例如：
12345# 运行前chr1 10000 10001chr1 10001 10002# 运行后chr1 10000 10002
但在最近遇到的分析场景中，我仅想合并重叠的点并对它们的值求平均，并不想将它们与相邻的点合并起来，例如：
123chr1 10000 10001 0.1chr1 10000 10001 0.2chr1 10001 10002 0.3
这里我仅想在 chr1:10000-10001 上进行合并且求平均，但此时使用 bedtools merge -c 4 -o mean 会输出以下结果：
1chr1 10000 10002 0.2
这并不符合我的预期工作情况，在查阅参数说明后，发现了以下关键参数：
-d：Ma ...
基于 US-align 计算两个蛋白结构之间的 RMSD 和 TM-score
在蛋白质结构研究里，比较不同蛋白质结构之间的相似性是一个核心任务，尤其是在结构预测、同源建模和结构功能关系研究中。以下是一些常用到的衡量蛋白结构间相似性的指标：
RMSD（均方根偏差）：该值是衡量两种蛋白质三维结构之间原子位置差异的标准度量。较小的 RMSD 值意味着两个结构相似度较高，RMSD 值对大分子在全局对齐时的大小和形状变化非常敏感。
TM-score：相较于 RMSD，TM-score 更加关注结构的全局相似性，且对于蛋白质的尺寸和形状变化具有更高的鲁棒性。TM-score 的值介于 0 到 1 之间，值越接近 1，表示两个结构越相似。
本文将介绍如何使用 US-align 对两个蛋白结构进行 align（对齐）并计算其 RMSD 及 TM-score。
US-alignUS-align github page：
https://github.com/pylelab/USalign
该软件的研究团队也同样是 TM-align 的开发者。相比其他的计算工具，US-align 有以下优点：
支持比对除蛋白质以外的生物大分子及其复合物（例如 DNA / RNA 等）。
通过 ...
使用 ESM 计算蛋白序列的突变影响分数
ESM（Evolutionary Scale Modeling） 是由 Facebook AI Research 开发的蛋白质语言模型。其基于大规模的蛋白质序列数据进行训练，能够捕捉蛋白质序列中的复杂模式和进化信息，从而在多个生物学下游应用中表现出色。
预训练好的 ESM 可用于多个生物学下游应用，本文将就其突变影响分析方面的功能进行介绍。
关于 ESMESM 预测突变影响的方式是无监督的，也就是说 ESM 的训练过程中并未使用到与突变影响相关的标签进行有监督学习。其评估每个氨基酸错义突变影响程度的方式是掩蔽该位置上的氨基酸并计算两种氨基酸（参考氨基酸与突变氨基酸）在该位置出现的对数似然比：
ESM 的预训练阶段使用了类似于 BERT 的掩码任务，因此 ESM 可以给出每个位点里各类氨基酸出现的概率，而参考氨基酸（Ref，即未突变时的氨基酸）和突变氨基酸（Alt）的概率差异即可作为一种 effect score。可以简单理解为，突变氨基酸出现的概率比参考氨基酸低越多，说明这类突变越有害。这种无监督的评估方式也被用于许多其他的深度学习工具，例如 EVE 和 GPN 等。
Refe ...
基于 JCVI 共线性分析确定一对一直系同源基因
前段时间因为课题需求，需要获取两个物种间的所有一对一直系同源基因对（One-to-One ortholog, 为避免文章内容冗余后面统称 OTO）。由于两个目标物种在 Ensembl 上都有，所以可以直接通过 BioMart 下载其基因的同源关系并提取出来。
但经过筛选和过滤后，发现 BioMart 中能确定的 OTO 数量极低，不太能满足分析需求。考虑到我所分析的目标基因类型在不同物种间具有高顺序保守性和序列保守性，因此想看看能不能通过 JCVI 的共线性结果来获取更多的 OTO 作为补充，发现效果还挺不错，固有此文。
JCVI 介绍JCVI 是一个用于比较基因组学分析的多功能工具包，共线性分析只是它其中的一个子功能。
JCVI 的安装，最简单的方式是通过 pip：
1pip install jcvi
也可通过 conda 或 mamba 安装：
1conda install bioconda::jcvi
此外也需要安装用于比对的 LAST：
1conda install bioconda::last
关于 JCVI 的运行，更多细节可见：JCVI tutorial
请注意， ...
获取 vcf 文件中每个个体携带的突变数量
最近的工作涉及到处理大量种群数据，其中有个分析需要统计每个个体携带的突变数量，在网上搜寻了诸多方法后确定了一些比较高效的策略，固有此文记录。
本文的部分流程参考源：
BioStars "how to count variants par sample per chromosome in a vcf file?"
提取每个个体携带的突变数量该流程主要通过 plink2 实现：
12plink2 --vcf [vcf file] --out [plink output1]plink2 --pfile [plink output1] --sample-counts --out [plink output2]
[vcf file]：此处填入 vcf 文件的路径。
[plink output1]：此处填入 plink2 输出的文件前缀。
[plink output2]：此处填入 plink2 输出的文件前缀。
运行结束后，可以得到以 .scount 结尾的文件，其中各列的含义如下：
#IID：标识每个样本的编号。
HOM_REF_CT：纯合参考基因型计数，个体纯合参考基因型的位点数量 ...
vscode 在远程服务器上通过 jupyter notebook 使用 R
最近一段时间的关键分析要用到 R，但是办公电脑的机械硬盘前段时间出了问题，现在用 Rstudio 会非常卡顿，加上相关分析对计算资源的要求也比较高，就想着是否有什么办法能在服务器上也弄一个类似于 Rstudio 的可视化界面出来作为辅助。
不过看了一圈发现 Rstudio 的 server 版本似乎局限性不少，而且还要用到 root 权限。考虑到平时运行 Python 进行分析时大多数时候都在 jupyter notebook 上进行，就想着是否能套用过来，于是就有了这篇文章。
安装该文章的具体实现在 Linux + mamba 环境下完成，如果你没有 conda/mamba 环境，可以参考网上资料进行安装。
首先确认你的环境中有 R，如果没有可以通过以下命令安装：
12mamba install conda-forge::r-baseR # 确认可以进入 R 终端，使用 q() 退出
通过 mamba 安装 IRkernel 和 jupyter（使用 conda 也行，本文章以 mamba 为准）：
12mamba install r::r-irkernelmamba instal ...
簇突变的类型简介及识别方法（Clustered mutation classification）
基因组中，某些特殊的突变过程可能会导致一小段区域内连续的突变发生。这些突变在基因组中呈现成簇分布，因此被称作 clustered mutation（本文将称其为簇突变）。研究这些突变的分布及组成有助于我们揭示导致其发生的内源性和外源性过程。该文章将主要讲述簇突变的类别并列举一些会导致特定类别的生物学原因，同时介绍该领域中的某些生物信息学工具等。
簇突变类型簇突变类型主要被分作六类，其中五类针对碱基水平的变化，一类针对插入及缺失（InDel）：
DSB（Doublet Base Substitutions）：双碱基替换，表示这两个突变发生在相邻的碱基上。某些外源性过程例如紫外线损伤会导致双碱基替换的发生（CC>TT），此外 DNA 修复缺陷和聚合酶功能突变也会导致 DSB 发生。
MBS (Multiple Base Substitutions)：多碱基替换，表示在很短的序列范围内发生多个碱基突变且这些突变彼此相邻。由于该簇突变的出现数量很有限，因此尚未得到全面研究。
Omikli：源自希腊语，意为 雾 或 薄雾，也被称作 diffuse hypermutation（弥漫性超 ...
rsync 中排除文件的规则详解（rsync --exclude）
前段时间在备份服务器时，由于备份策略的变化，对 rsync 的 exclude patterns 进行了一些更加深入的探索。
这篇文章将介绍如何设计 rsync 的 exclude patterns 以排除特定目录或文件。
命令介绍rsync 是 Linux 和 Unix 中常用的文件复制和同步命令，相比 cp 命令而言更加高效和完善，其优点包括但不限于：
增量复制：rsync 只会复制那些自上次同步后更改的文件部分，极大地减少了数据传输量。例如 a 文件已经同步过一次并且没有发生变化，则 rsync 会跳过。同样，rsync 也会进行完整性检查以确定传输的数据无误。
保留属性：rsync 可以保持文件的权限、时间戳、软硬链接、所有权等属性。
灵活的过滤规则：支持使用包括和排除规则来选择同步的文件或目录。
排除规则rsync 通过参数 --exclude 或 --exclude-from 指定排除的文件模式。其中 --exclude-from 用于指定包含排除模式的文件（包含多个排除规则）。
排除指定文件输入文件名称（如 123.txt）即可排除所有命名为 123.txt 的文 ...
Snakemake pipeline 搭建的进阶教程
关于 SnakemakeSnakemake 是一个工作流管理系统，可用于创建易于迁移和复现的数据分析流程（pipeline），因此被广泛应用于生物学分析中。虽然 Snakemake 算是一个独立的编程语言，但其本质是基于 Python 搭建的，因此如果你对 Python 有一定了解，那么上手 Snakemake 就会更容易很多。并且 Python 技能的掌握对于某些情况下 Pipeline 下游分析部分的拓展是有必要的。
本文将结合多个近些年生物论文中公开的 Snakemake Pipeline，介绍如何搭建一个具有一定功能的 Pipeline。
Snakemake 安装：https://snakemake.readthedocs.io/en/stable/getting_started/installation.html
Snakemake Pipeline 示例教程：https://snakemake.readthedocs.io/en/stable/tutorial/basics.html
有利于提升阅读体验的一些条件：
对于 Python 编程已经具备一定的了解，包括文件 ...
对 bedgraph 文件进行 lowess 平滑操作
在进行一些特定的功能基因组学数据分析时，我们可能需要对 bedgraph 文件中每个 bin 的值进行一定的平滑操作，以降低随机噪声的影响并提供更好的可视化效果。例如：
Repli-seq / BrdU-seq 中量化得到的 Replication Timing
OK-seq / Pu-seq 中量化得到的 Replication Fork Directionality
以下是一个用于对 bedgraph 进行 lowess 平滑操作的 python script：
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960import pandas as pdimport statsmodels.api as smimport argparsedef get_args(): parser = argparse.ArgumentParser(description='Perform lowess smooth on ...
获取基因组上 intron 区域及对应链信息的方法
流程前置条件：
在环境变量中可调用的 bedtools
本文参考：
Get intronic and intergenic sequences based on gff file from Biostars
https://www.biostars.org/p/112251/
本文可满足的需求：
得到特定区域的 bed 文件（例如 exon / intron 等）。
在得到区域信息的同时进行链信息的区分。
下载基因组注释文件（gtf）以人类最新版本的 gencode 注释为例：
1wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_46/gencode.v46.basic.annotation.gtf.gz
也可以选择已有的 gtf 文件进行。
得到 transcript 和 exon 区域信息12345678910111213141516GTF="/path/to/gencode.v46.basic.annotation.gtf.gz" # 改为自己实际的注释文件路径，建议使用压缩 ...
Ensembl VEP plugins 的使用方法（Alphamissense、dbNSFP 等）
前言以下是相关的官网链接，如果想要得到除本文内容以外更全面深入的了解，建议跳转相关链接进行查阅：
VEP plugins documentation:
https://grch37.ensembl.org/info/docs/tools/vep/script/vep_plugins.html
VEP install documentation:
https://grch37.ensembl.org/info/docs/tools/vep/script/vep_download.html
正文安装 VEP此处强烈建议使用 Docker 或者 Singularity 直接拉取 VEP 镜像进行分析，因为在 VEP image 中已经内置好所有的 plugin，无需自己手动下载。
考虑到 Docker 对权限的要求更高，因此这里以普适性更强的 Singularity 为例（Docker 用户有需求可通过文章开头链接前往官网进行参考）：
12# 拉取 VEP imagesingularity pull --name vep.sif docker://ensemblorg/ensembl- ...