提取 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 进行过滤， ...
bismark 分析 Bisulfite-Seq 数据的流程示例
前言以下是相关的官网链接和参考资料，如果想要得到除本文内容以外更全面深入的了解，建议直接跳转相关链接进行查阅：
Bismark github:
https://github.com/FelixKrueger/Bismark
Bismark documentation:
https://felixkrueger.github.io/Bismark/
本文比对及去重等流程均参考自 Bismark documentation
此外，以下内容在本文中不会涉及到，有需要请自行探索：
数据的质控
除 Bismark 以外其他软件的安装
通过以下流程可以得到：
基因组上 CpG 位点的甲基化信息
请阅读完全文后再进行相关分析，建议先使用 bismark 的示例数据观察相关命令是否存在问题。
正文下载最新版本的 bismark（至 2024/04/19）：
1$ wget https://github.com/FelixKrueger/Bismark/archive/refs/tags/v0.24.2.tar.gz
也可前往 github release 下载指定版本或可能已存在的更新 ...
通过 bowtie2 + macs3 分析 ChIP-Seq 数据
前言以下是相关的官网链接和参考资料，如果想要得到除本文内容以外更全面深入的了解，建议直接跳转相关链接进行查阅：
MACS3 github:
https://github.com/macs3-project/MACS
MACS3 documentation:
https://macs3-project.github.io/MACS/
本文比对及去重等流程参考自以下资料：
Alignment and filtering from Harvard Chan Bioinformatics Core (HBC)
Histone ChIP-seq pipeline from ENCODE project (Version 2.0 with bowtie2)
本文流程根据以下资料进行过优化：
sambamba -F “not duplicate” processed bam still have duplicated (sambamba issue#477) (Update 2024/4/24)
Encode blacklist (Update 2024/4/24)
此外，以下内容在本 ...
远程服务器无 root 权限时 LightGBM GPU 版的安装方法
正文首先感谢 LightGBM 团队的帮助：
https://github.com/microsoft/LightGBM/issues/6399
声明本文适用于在远程服务器（Linux）上安装 LightGBM 的情况。
不同的编译器版本在遇到报错时的解决方案可能略有不同，以下是我在编译时使用的相关版本：
12-- The C compiler identification is GNU 8.5.0-- The CXX compiler identification is GNU 4.8.5
本文默认已在服务器上配置有 CUDA，如果 CUDA 另有路径请通过 -DOpenCL_LIBRARY 和 -DOpenCL_INCLUDE_DIR 指定。
编译下载 LightGBM，并进行 CUDA 版本的编译：
123456git clone --recursive https://github.com/microsoft/LightGBMcd LightGBMmkdir buildcd buildcmake -DUSE_CUDA=1 .. # 指定 CUDA versionmake -j ...
全基因组 CpG (or nonCpG) 位点 bed 文件的获取方式
CpG 位点是基因组上的一种特殊序列，它在 DNA 的结构和功能中起着重要作用。其组成和甲基化影响着基因组上各区域的调控，是表观遗传学领域中重要的研究目标。
本文中的脚本主要用于提取基因组中 CpG 位点的位置，以方便一些着重于 CpG 位点的下游分析进行。此外在部分场景下，我们可能希望排除 CpG 位点进行分析（因为其独特的性质例如高突变率等），因此该脚本也可产生 AT&CG(nonCpG) 位点信息文件。
需要准备的文件和条件有：
基因组序列文件
Biopython 库（pip install biopython）
可能需要准备的前置条件：
已被放置于 $PATH 环境变量中的 bedtools
脚本内容：
1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677from Bio import SeqIOimport osimport ...
Clash for windows 中让某些域名绕过代理
众所周知，Clash for windows 的作者已经删库跑路许久，而在使用这个软件到现在我也终于忍不住想要解决某些域名绕过代理的问题。我已经受够了每次进 web of science 查文献亦或是在 Science 官网上看文章都要关掉代理才能进去的痛楚了。
在搜索网上的解决方案时，部分文章说可以在 Settings>System Proxy>Bypass Domain/IPNet 中设置相应域名解决，但是我发现我新加域名进入后我的梯子就会直接整个没用。当然这发生在大约半年前的时候，那时的我选择了将这件事暂时搁置，直到现在。
这次再次寻找解决方法，终于找到个适用于我的，趁这个机会写在此处，希望能帮到其他有需求的朋友，当然也有可能某天这篇文章会因为不可抗力消失。
解决方案首先点开 Profiles，选择自己的梯子右键，点开 Parsers 并选择 Edit Parsers。
将里面的内容处理成以下格式：
12345parsers: # array - url: 梯子 url yaml: prepend-rules: - DOM ...
通俗易懂的贝叶斯深度神经网络介绍及代码实现
最近在课题组里的一些任务涉及到了把模型的架构改成贝叶斯层，让模型输出的结果能够体现不确定性。恶补相关知识的同时正好也能记录下笔记发博客，一举两得。
这篇文章只做概念介绍，不做深入剖析，追求原理的话敬请走下文中相关链接前往原论文。
贝叶斯神经网络是什么或许你在很多其他的地方听过贝叶斯这个词汇。以下是一些生物学中的例子：
系统发育学中的贝叶斯建树（Mrbayes、Phylobayes）。
群体遗传学（种群遗传学）中分析种群结构、估计种群历史等（Structure、ABC）。
在某些基因组比对的工具里也有用到贝叶斯方法（bwa）。
贝叶斯方法的大致原理是结合先验知识和实验数据来估计参数或测试假设，这个过程涉及到随机采样，因此其他方法不同，应用贝叶斯法得到的结果往往不会是完全一致的，也就是说，可能每次跑出来的结果之间会相似，但不会完全相同。
这里就提及到了贝叶斯方法的一个宝贵性质 —— 不确定性。
如果有模型搭建经验的话，你应该知道，对于传统的神经网络架构而言，在你训练好模型后，如果给它相同的输入，那么一般而言它都会产生一致的输出（因为模型内部的权重是确定的）。但是有些时候，我们可能想知道 ...
dN/dS 的计算及其核心理念
dN/dS (ω) 是一个在进化生物学领域常出现的词，也是我最近分析的时候着重关注的一个对象。这几天翻了一些文献和书籍，对它的计算和设计理念有了一个全新的了解，所以有了这篇博客。
在阅读这篇文章前，你可能需要一定的相关知识基础，例如，你应该已经知道以下概念：
密码子、简并位点、同义突变、非同义突变。
中性理论、正选择、纯化选择。
dN/dS 是什么这里的 S 即 Synonymous，N 即 Nonsynonymous。
dN 指的是非同义替换率，dS 指的是同义替换率。但我认为另一个知乎博主文章里的说法会更加直观一些：
dN: 每个非同义位点上发生的非同义突变的数量。
dS: 每个同义位点上发生的同义突变的数量。
计算原理首先，简单来说，dN 和 dS 的计算涉及到四个变量，分别为同义突变数量（S）、非同义突变数量（N）以及同义位点数量（Ssite）和非同义位点数量（Nsite）。
前两个变量 N 和 S 非常好理解，所以这里需要重点关注的是后两个变量 Nsite 和 Ssite，即每种类型的位点数量是什么？它们是如何得到的？但在这之前，我们先简单回顾一下一些常见的突变类 ...
slurm 作业数组的正确使用姿势
最近一次想要批量运行任务前试着运行了一个个例，发现占用的内存意外的比较多，所以一口气把所有提交到服务器上显然不太妥当。由于之前在服务器上用 slurm 仅限于 srun 而从来没有尝试过 sbatch，所以就对这方面的知识进行了一番恶补，正好也能补充十二月份的博客内容。
如何确定内存是否充足对于当前节点，使用 top 和 htop 就能得到相关的信息，但对于其他节点，如果没有登录权限，那么这种实时监控的手段就没法派上用场了。
这时候，可以使用 free -h 命令将节点在该时间点的内存使用情况信息打印出来：
1234$ srun -w nodexx free -h total used free shared buff/cache availableMem: 487G 144G 2.7G 162M 340G 340GSwap: 49G 685M 49G
重点关注 availa ...
DNA 语言模型 GPN 在突变效应预测中的应用及表现
众所周知，语言模型已经在诸多领域得到了广泛应用。受深度学习狂潮影响深远的生物学领域里，出色的语言模型应用不断涌出。以 DeepMind 为例，其所开发的 Alphafold、Alphamissense 等都是基于语言模型实现的。
这篇文章将分享一个最近发表的 DNA 语言模型（GPN，Genomic Pre-trained Network）和它的变体 GPN-MSA，该模型与之前的 DNABERT 等 DNA 语言模型不同，其使用单个碱基作为 token 进行训练，并在无监督的突变效应预测中取得了良好的成效。
相关论文链接见文章结尾，本文章将重点介绍两个模型的设计和它们之间的区别，并写出一些很值得学习的地方。
GPNGPN 是一种自监督训练的 DNA 语言模型，与 GPN-MSA 不同的地方在于，它可以仅依靠未经过比对的数个基因组序列而不使用多序列比对（MSA）进行训练，此外它使用的为单纯的 DNA 序列而未使用到任何功能基因组信息，因此可以说 GPN 的泛用性是非常强大的，因为对于很多类群中的生物来说，获得跨物种的 MSA 和全面的功能基因组数据并不是容易的事情。
GPN 的模型设计 ...
深度学习 loss 追踪器（用于追踪 loss 变化）
最近设计了一个深度神经网络，虽然使用 logging 模块记录 Train loss 和 Valid loss 也算方便，但作图观察 loss 的变动趋势还是要更加直观的，考虑到网上貌似没有这方面的方法，这里放一下自己的记录方式。
1234567891011121314151617181920212223242526272829303132import matplotlib.pyplot as pltfrom matplotlib.ticker import MaxNLocatorimport timeclass loss_tracer(): def __init__(self, outpath): self.train_loss_list = [] self.test_loss_list = [] self.outpath = outpath timestamp = time.time() self.starttime = time.strftime('%Y_%m_%d_%H_%M', time.localtime(timestamp)) def __call__(sel ...
Python 循环太费时？让多进程解决这个问题（concurrent）
前言这篇文章主要分享如何使用 Python 的 concurrent 库进行多进程运行以加速分析。在约一年前我在跑多个物种的 ABBA-BABA test 时遇到过分析缓慢的问题，因为涉及到的物种有九十多种，所以在经过筛选以后四物种组合的数量依然高达数百万，但由于那个分析并不是很要紧所以也就没有多管（最后跑了一个多月）。
现在的状况与那时有些许不同，最近跑的分析大多涉及到大量计算，由于比较急迫地想看到最终结果，我不得不学习多进程方法来进行加速，很多时候进步都需要一些些小小的契机，now it comes。
由于 Python 并不是从底层框架开始让人写代码，所以它在拥有友好入门门槛的同时也牺牲了一定的性能，特别是在处理计算密集型任务时，Python 比起 C 来说就显得有些捉襟见肘。
当然我肯定是不愿意写 C 的代码的，一是因为我只学了一些皮毛，二是因为相比之下 Python 实在太方便，俗话说性能不够进程来凑，既然我跑的不够快，那我就多跑一些，以量胜质。
在使用多进程前，建议先检查一下自己的代码是否可以进行优化，有时候换成合适的数据操作方式就可以让分析的速度有质的飞跃，使用多进程应 ...
Python seaborn 绘制复杂多变量图（hexbin + 拟合直线）
当初作图时想往 sns.pairplot 里加些内容，搜索了解到 sns.PairGrid 可以更灵活地绘制多因素（多变量）关系图，遂在网上找相关教程，最后发现 CSDN 上的内容竟然都要付费，无奈转向 ChatGPT 求助，最后在 C 老师帮助下成功完成了自己的想法。
以上既是这篇文章的来源也是动机，这篇文章就来分享一下如何使用 sns.PairGrid 绘制更加复杂的多变量图，本文示例将是 hexbin + 拟合直线，但只要知道如何将图架构起来那么绘制其他类型的图也不会有问题。
sns.pairplot 的详尽用法：
https://seaborn.pydata.org/generated/seaborn.pairplot.html#seaborn-pairplot
一般的多变量图绘制：
12345678910111213141516171819import numpy as npimport pandas as pdimport matplotlib.pyplot as pltimport seaborn as sns# 模拟三个变量np.random.seed(114514)a ...
关于人类基因组注释中起始和终止位点相同的现象
前言已经有一个月没有更新文章，之前在知乎看到过一段话，大体是说写博客最重要的不是怎么把博客做的好看，而是如何持之以恒地产出，这句话让我感触颇深，正值国庆期间放假，就记录一下这段时间在分析数据时候发现的东西，以后万一有其他朋友发现并疑惑也好快速找到答案。
基因组注释中长度为 1 的现象解释在使用基因组注释（版本：https://www.gencodegenes.org/human/release_19.html ）时，我发现了两个问题：
有些 CDS 的长度为 1。
CDS 的长度不一定为 3 的倍数。
由于我是先发现的第一个问题，此后才发现的第二个问题，因此顺序如上，首先看一看第一个问题：
12345678910111213# First, download the GENCODE v19 annotation file$ wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz# Next, extract and count t ...