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- ...
fasterq-dump 下载 SRA 文件时报错的解决方法（err cmn_iter）
正文这几天下载 SRA，遇到的错误有：
fasterq-dump.3.1.0 err: cmn_iter.c cmn_read_uint8_array
fasterq-dump.3.1.0 err: cmn_iter.c cmn_read_String
以前使用 fasterq-dump 时只是偶然出现这些问题，重新下载也都能解决，但最近有些 SRA 一直下载失败，经过调查发现这种错误在大文件下载时出现异常频繁，且多次下载并没法有效解决问题，因此需要一个替代的方法。
解决方案Prefetch由于 fasterq-dump 直接通过 HTTP 下载得到 fastq 文件，该过程很可能由于一些问题中断从而导致下载失败。因此可以通过更稳定的 prefetch 先得到 sra 文件，再通过 fasterq-dump 提取 fastq 文件。
fasterq-dump fetches SRR on the fly via HTTP and there could be fatal errors during the transfer.prefetch eliminates transf ...
提取 bwa 比对中唯一比对 reads（uniquely mapped reads）的方法
本文由 Juse 基于以下资料进行撰写：
How to extract uniquely mapped reads from bam/sam produced by bwa-mem? from SEQanswers
Obtaining uniquely mapped reads from BWA mem alignment (filtering by q score does not seem to do the trick in my case) from Biostars
Obtaining uniquely mapped reads from BWA mem alignment from StackExchange
在此感谢 community 中各位大佬的无私分享。
正文在某些文章中，有时我们会看到作者提及到在后续分析中，ta 只考虑了 uniquely mapping reads（后文统称唯一比对）。从直觉上出发，使用唯一比对的思路不难理解 —— 这些比对更加精确，并使后续的分析具有更好的解释性。
目前来说，得到唯一比对的最主流方法是直接通过 MAPQ 进行过滤，但这并不 ...