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 进行过滤，但这并不 ...
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 ...