其他

在iPad上阅读文献和图书

长久以来，我用iPad读文献的操作一直是下面这几步

  1. 用Zotero从网上抓取文献
  2. 找出我的iPad连接线
  3. 将iPad连接到电脑
  4. 找到Zotero文献的所在目录
  5. 打开iTunes的共享
  6. 把PDF拖到iTnues的共享中
  7. 打开iPad用GoodReader打开文献

每次插连接线都感觉自己太过原始。后来有了有了一部iPhoen手机，于是就变成了

  1. 用Zotero从网上抓取文献
  2. 打开电脑微信，把PDF发给自己
  3. 打开iPhone的微信
  4. 用iPhone的AirDrop功能将PDF共享给我的iPad
  5. 打开iPad用GoodReader打开文献
  6. 清理电脑微信和手机微信上的文献

手动清理是在是蠢。好在后来有了坚果云，就变成了

  1. 用Zotero从网上抓取文献
  2. 把文献保存到坚果云
  3. 在iPad的GoodReader上同步文献

但是GoodReader这个软件只能作为一个阅读工具,它的PDF编辑功能太弱, 也不适合作为一个笔记工具

现在知道了PaperShip和MarginNote 3，于是就是

  1. 用Zotero抓取文献
  2. 用PaperShip同步文献，阅读文献
  3. 假如文献信息量很大，用MarginNote 3做笔记
  4. MarginNote 3做的笔记可以导出至印象笔记

不过PaperShip这个软件存在一些问题，会出现文件无法同步的问题，只能期待有更好的App了。感谢坚果云，目前PaperShip能够顺利进行同步了。

只要打开坚果云官网，登录坚果云账号，在zotero文件夹新建空白的lastsync.txt文件，务必注意的是，必须使用坚果云自带的新建文件工具来新建lastsync.txt文件，不能通过手动上传的方式。

新建一个文件

接着，重新在Zotero客户端或者Papership验证服务器即可

验证成功

如果我想看的是一本书，那么就是把书放到坚果云中，用MarginNote 3打开坚果云里的书做笔记。

利用ZotFile进行文献同步批注

还有一个解决方案，是ZotFile这款插件。我们只需要新建一个坚果云同步文件夹，Tablet，并在Zotfile的插件中加入该选项。

同步文件夹

此时在Zotero处就会多出两个层次，存放送去批注的文件，以及批注完成的文件。

Tablet Files

我们将要阅读文献右击，在Manage Attachment 中选择 Send to Tablet

送至平板

因为我在GoodReader上设置坚果云中的Tablet作为同步服务器，因此就可以在GoodReader对文献进行批注。批注完成后，在Tablet Files(modified)中选择批注完成的文件，右击，在Manage Attachment 中选择 Get From Tablet

取回修改后内容

感觉这个方法也很不错，毕竟你回一次寝室看的文献也不会有那么多。

参考资料

其他

使用坚果云作为Zotero的存储

Zotero的存储2GB就需要20美刀，差不多的价格(199.90)能买一年的专业版坚果云，容量是54GB.

更重要的一点是用坚果云作为Zotero的同步服务器很方便。

第一步: 在坚果云上设置应用密码

1应用密码

第二步: 打开Zotero的首选项（preference），有可能是在工具栏中，也有可能是编辑栏中。

首选项

第三步: 在Zotero上关联坚果云的网盘，需要注意的是此处的密码不是坚果云的登录密码，而是第三方应用授权时的密码。

配置同步

后面同步就行了。

参考资料

官方教程: http://help.jianguoyun.com/?p=3168

Linux

服务器硬盘管理笔记

服务器硬盘分区和挂载

文件系统简介

以ext4文件系统为例，设计的时候分为4个部分

  • 超级块: 最开头部分，记录整个文件系统/分区中的文件数，df根据该信息进行统计磁盘占用
  • 超级块副本: 对超级块进行备份，存在多个备份
  • i节点(inode): 记录每个文件信息，名称/大小/编号/权限。 用ls -i获取i节点信息
  • 数据块(datablock): 事实上的数据存储点，i节点以链接式结构记录数据块的信息，因此找到i节点，就能读取对应的数据。

由于ls -l获取的是i节点记录的数据使用的数据块个数，而du则是通过i节点获取实际大小, 所以ls -l和du显示的数据大小不同。

RAID

RAID全称是Redundant Array of Independent Disks，也就是磁盘阵列，通过整合多块硬盘从而提升服务器数据的安全性，以及提高数据处理时的I/O性能。

RAID目前常用的是RAID5, 至少需要3块硬盘，其中一块硬盘用于奇偶校验，保证数据安全，其余硬盘同时读写，提高性能。此外，你还需要知道最原始的是RAID0，同时将数据读写到所有硬盘里，速度就变成了原来的N倍。RAID1至少需要两块盘，其中一块硬盘是另外硬盘的镜像。它不提高读写效率，只提高了数据安全性。RAID10是RAID0和RAID1的组合。

目前的服务器都配备了硬件RAID卡，因此在为服务器增加或更换硬盘时，需要格外注意，

  • 插入的硬盘必须先做RAID才能被服务器识别
  • RAID5中的所有硬盘应当被视为一个硬盘，不能只更换其中一部分

硬盘分区

fdisk只能对不多于2TB的硬盘进行分区

1
fdisk /dev/sdb

假如你的硬盘大于2TB，那么会输出如下信息

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
Welcome to fdisk (util-linux 2.23.2).

Changes will remain in memory only, until you decide to write them.
Be careful before using the write command.

Device does not contain a recognized partition table
Building a new DOS disklabel with disk identifier 0x405fb398.

WARNING: The size of this disk is 10.0 TB (10000294477824 bytes).
DOS partition table format can not be used on drives for volumes
larger than (2199023255040 bytes) for 512-byte sectors. Use parted(1) and GUID
partition table format (GPT).


The device presents a logical sector size that is smaller than
the physical sector size. Aligning to a physical sector (or optimal
I/O) size boundary is recommended, or performance may be impacted.

Command (m for help):

提示信息中的警告中，就建议”Use parted(1) and GUID partition table format (GPT).”

因此，对于大于2TB的硬盘就需要用parted进行分区

1
parted /dev/sdb

输出信息如下

1
2
3
GNU Parted 3.1
Using /dev/sdb
Welcome to GNU Parted! Type 'help' to view a list of commands.

创建新的GPT标签，例如

1
mklabel gpt

设置单位

1
unit TB

创建分区, 比如我将原来的10T分成2TB和8TB

1
2
3
# mkpart PART-TYPE [FS-TYPE] START END
mkpart primary  ext4 0.00TB 2.00TB
mkpart primary  ext4 2.00TB 10.00TB

查看分区表

1
print

输出如下

1
2
3
4
5
6
7
8
9
10
11

Model: AVAGO MR9361-8i (scsi)
Disk /dev/sdb: 10.0TB
Sector size (logical/physical): 512B/4096B
Partition Table: gpt
Disk Flags: 

Number  Start   End     Size    File system  Name     Flags
 1      1049kB  2000GB  2000GB               primary
 2      2000GB  10.0TB  8000GB               primary

退出

1
quit

此时会提示”Information: You may need to update /etc/fstab.” /etc/fstab用于设置开机硬盘自动挂载。如果硬盘被拔走了，而/etc/fstab没有修改，那么会就提示进行修复模式。

硬盘格式化

在挂载硬盘之前，需要先对磁盘进行格式化。使用的命令为mkfs, 使用-t指定文件系统，或者用mkfs.xxx，其中xxx就是对应的文件系统。文件系统有如下几类

  • 传统文件系统: ext2, minix, msdos, fat, vfat
  • 日志文件系统: ext3, xfs, ext4
  • 网络文件系统: nfs

目前最流行的是ext4和xfs，足够稳定。其中xfs是CentOS7之后的默认文件系统。

1
2
mkfs.ext4 /dev/sdb1
mkfs.ext4 /dev/sdb2

硬盘挂载

之后用mount进行硬盘挂载，分别两种情况考虑

一种是新建一个文件路径，进行挂载。

1
2
mkdir data2
mount /dev/sdb2 /data2/

另一种是挂载一个已有目录，比如说临时文件目录/tmp挂载到新的设备中。

第一步: 新建一个挂载点，将原有数据移动到该目录下

1
2
3
mkdir /storage
mount /dev/sdb1  /storage
cp -pdr /tmp/* /storage

第二步: 删除原来的/tmp下内容

1
rm -rf /tmp/*

第三步: 重新挂载

1
2
umount /dev/sdb1
mount /dev/sdb1 /tmp

和mount相关的文件如下

  • /etc/fstab: 文件系统表， 开机的时候会根据里面的记录进行硬盘挂载
  • /etc/mtab: 记录着已经挂载的文件系统
  • /etc/mtab~: 锁定文件
  • /etc/mtab.tmp: 临时文件
  • /etc/filesystems: 系统支持的文件系统

此外mount在挂载的时候还可以设置文件系统参数，例如是否支持磁盘配额，对应-o参数

总结

第零步: 检查服务器是否具备RAID阵列卡，如果有，则需要先为硬盘做RAID。

第一步: 使用fdisk -l检查硬盘是否能被系统检测到

第二步(可选): 假如需要硬盘分区，则用fdisk/gdisk/parted对硬盘划分磁盘

第三步: 使用mkfs进行磁盘格式化，有如下几种可选，

第四步: 用mkdir新建一个目录，然后用mount将格式化的硬盘挂载到指定目录下。卸载硬盘，则是umout

第五步: 修改/etc/fstab 将硬盘在重启的时候自动挂载。注意: 如果硬盘不在了，则需要将对应行注释掉，否则会进入到emergency模式。

生信软件工具箱

使用OrthoFinder进行基因家族分析

谈论到直系同源基因分析的时候，大部分教程都是介绍OrthoMCL，这是2003年发表的一个工具，目前的引用次数已经达到了3000多，但这个软件似乎在2013年之后就不在更新，而且安装时还需要用到MySQL（GitHub上有人尝试从MySQL转到sqlite）。

而OrthoFinder则是2015年出现的软件，目前已有400多引用。该软件持续更新，安装更加友好，因此我决定使用它来做直系同源基因的相关分析。

OrthoFinder能做什么？

OrthoFinder: solving fundamental biases in whole genome comparisons dramatically improves orthogroup inference accuracy提到，它的优点就是比其他的直系同源基因组的推断软件准确，并且速度还快。

此外他还能分析所提供物种的系统发育树，将基因树中的基因重复事件映射到物种树的分支上，还提供了一些比较基因组学中的统计结果。

OrthoFinder的分析过程

OrthoFinder的分析过程分为如下几步:

  1. BLAST all-vs-all搜索。使用BLASTP以evalue=10e-3进行搜索，寻找潜在的同源基因。(除了BLAST, 还可以选择DIAMOND和MMSeq2)
  2. 基于基因长度和系统发育距离对BLAST bit得分进行标准化。
  3. 使用RBNHs确定同源组序列性相似度的阈值
  4. 构建直系同源组图(orthogroup graph)，用作MCL的输入
  5. 使用MCL对基因进行聚类，划分直系同源组

分析流程1
OrthoFinder2在OrthoFinder的基础上增加了物种系统发育树的构建，流程如下

  1. 为每个直系同源组构建基因系统发育树
  2. 使用STAG算法从无根基因树上构建无根物种树
  3. 使用STRIDE算法构建有根物种树
  4. 有根物种树进一步辅助构建有根基因树

基于Duplication-Loss-Coalescent 模型，有根基因树可以用来推断物种形成和基因复制事件，最后记录在统计信息中。

分析流程2

软件使用

在解压缩的OrthoFinder文件目录下(安装见最后)有一个 ExampleData, 里面就是用于测试的数据集。

1
2
3
orthofinder -f ExampleData -S mmseqs
# -f 指定文件夹
# -S 指定序列搜索程序，有blast, mmseqs, blast_gz, diamond可用

OrthoFinder的基本使用就是如此简单，而且最终效果也基本符合需求。

如果你想根据多序列联配(MSA)结果按照极大似然法构建系统发育树，那么你需要加上-M msa。这样结果会更加准确，但是代价就是运行时间会更久，这是因为OrthoFinder要做10,000 - 20,000个基因树的推断。

OrthoFinder默认用mafft进行多序列联配，用fasttree进行进化树推断。多序列联配软件还支持muscle, 进化树推断软件还支持iqtree, raxml-ng, raxml。例如参数可以设置为-M msa -A mafft -T raxml.

并行化参数: -t参数指定序列搜索时的线程数，-a指的是序列搜索后分析的CPU数。

软件细节

OrthoFinder提供了config.json可以调整不同软件的参数，如下是BLASTP。

BLASTP参数

OrthoFinder默认使用DendroBLAST发育树，也就是根据序列相似度推断进化关系。这是作者推荐的方法，在损失部分准确性的前提下提高了运算效率。当然你可以用-M msa从多序列比对的基础上进行基因树构建。如果你先用了默认的DendroBLAST，想测试下传统的MSA方法，那么也不需要重头运行，因为有一个-b参数可以在复用之前的比对结果。

在物种发育树的推断上，OrthoFinder使用STAG算法，利用所有进行构建系统发育树，而非单拷贝基因。此外当使用MSA方法进行系统发育树推断时，OrthoFinder为了保证有足够多的基因(大于100)用于分析，除了使用单拷贝基因外，还会挑选大部分是单拷贝基因的直系同源组。这些直系同源组的基因前后相连，用空缺字符表示缺失的基因，如果某一列存在多余50%的空缺字符，那么该列被剔除。最后基于用户指定的建树软件进行系统发育树构建。结果在”WorkingDirectory/SpeciesTree_unrooted.txt”

使用STRIDE算法从无根树中推断出有根树, 结果就是”SpeciesTree_rooted.txt”.

结果文件

运行结束后，会在ExampleData里多出一个文件夹，Results_Feb14, 其中Feb14是我运行的日期

直系同源组相关结果文件，将不同的直系同源基因进行分组

  • Orthogroups.csv：用制表符分隔的文件，每一行是直系同源基因组对应的基因。
  • Orthogroups.txt: 类似于Orthogroups.csv，只不过是OrhtoMCL的输出格式
  • Orthogroups_UnassignedGenes.csv: 格式同Orthogroups.csv，只不过是物种特异性的基因
  • Orthogroups.GeneCount.csv：格式同Orthogroups.csv, 只不过不再是基因名信息，而是以基因数。

直系同源相关文件，分析每个直系同源基因组里的直系同源基因之间关系，结果会在Orthologues_Feb14文件夹下，其中Feb14是日期

  • Gene_Trees: 每个直系同源基因基因组里的基因树
  • Recon_Gene_Trees：使用OrthoFinder duplication-loss coalescent 模型进行发育树推断
  • Potential_Rooted_Species_Trees: 可能的有根物种树
  • SpeciesTree_rooted.txt: 从所有包含STAG支持的直系同源组推断的STAG物种树
  • SpeciesTree_rooted_node_labels.txt: 同上，只不过多了一个标签信息，用于解释基因重复数据。

比较基因组学的相关结果文件：

  • Orthogroups_SpeciesOverlaps.csv： 不同物种间的同源基因的交集
  • SingleCopyOrthogroups.txt： 单基因拷贝组的编号
  • Statistics_Overall.csv：总体统计信息
  • Statistics_PerSpecies.csv：分物种统计信息

STAG是一种从所有基因推测物种树的算法，不同于使用单拷贝的直系同源基因进行进化树构建。

一些重要概念:

  • Species-specific orthogroup: 一个仅来源于一个物种的直系同源组
  • Single-copy orthogroup: 在直系同源组中，每个物种里面只有一个基因。我们会用单拷贝直系同源组里的基因推断物种树以及其他数据分析。
  • Unassigned gene: 无法和其他基因进行聚类的基因。
  • G50和O50，指的是当你直系同源组按照基因数从大到小进行排列，然后累加，当加入某个组后，累计基因数大于50%的总基因数，那么所需要的直系同源组的数目就是O50，该组的基因树就是G50.

Orthogroups, Orthologs 和 Paralogs 这三个概念推荐看图理解。

概念辨析

如何安装？

最快的方法

OrthoFinder可以通过conda安装，建议为它新建一个虚拟环境

1
conda create -n orthofinder orthofinder=2.2.7

如果你愿意折腾

你先得安装它的三个依赖工具: MCL, FastME, DIAMOND/MMseqs2/BLAST+

MCL有两种安装方式，最简单的就是用sudo apt-get install mcl, 但是对于大部分人可能没有root权限，因此这里用源代码编译。http://micans.org/mcl/

1
2
3
4
5
wget https://www.micans.org/mcl/src/mcl-latest.tar.gz
tar xf mcl-latest.tar.gz
cd mcl-14.137
 ./configure --prefix=~/opt/biosoft/mcl-14.137
make -j 20 && make install

之后是MMseqs2, 一个蛋白搜索和聚类工具集，相关文章发表在NBT, NC上。GitHub地址为https://github.com/soedinglab/MMseqs2

1
2
3
wget https://github.com/soedinglab/MMseqs2/releases/download/3-be8f6/MMseqs2-Linux-AVX2.tar.gz
tar xzf MMseqs2-Linux-AVX2.tar.gz
mv mmseqs2 ~/opt/biosoft/

最后安装FastME, 这是一个基于距离的系统发育树推断软件。在http://www.atgc-montpellier.fr/fastme/binaries.php下载，上传到服务器

下载

1
2
3
4
tar xf fastme-2.1.5.tar.gz
cd fastme-2.1.5
./configure --prefix=/opt/biosoft/fastme-2.1.5
make && make install

BLAST+可装可不装，推荐阅读这或许是我写的最全的BLAST教程

以上软件安装之后，都需要将其添加到环境变量中，才能被OrthoFinder调用。

之后在https://github.com/davidemms/OrthoFinder/releases 寻找最近的稳定版本下载到本地，例如OrthoFinder v2.2.7

1
2
tar xzf OrthoFinder-2.2.7.tar.gz
OrthoFinder-2.2.7/orthofinder -h
生信软件工具箱

使用OrthoMCL鉴定直系同源基因组

OrthoMCL是目前最常用的基因家族分析软件，从2013年发布2.0版本之后再也没有更新过，虽然它的安装过程复杂负责，但是依旧挡不住大家对他的喜好。

当然软件安装复杂是相对于之前，现在用Docker就可以轻松的安装和使用OrthoMCL。由于

1
docker pull jasonkwan/orthomcl_docker

由于该容器自带一个run_orthomcl.py, 只要你准备好输入数据和配置文件，就能够自动化进行分析.

其中run_orthomcl.py的参数有5个

  • -t/–table_path fasta_table的文件路径
  • -s/–start_stage: 开始阶段，1,2,3,4
  • -e/–end_stage: 结束阶段， 1,2,3,4
  • -p/–processors: 线程数
  • -c/–config_file: 配置文件

其中-s和-e的1-4分别对应为

  • 1: 使用orthomclAdjustFasta, orthomclFilterFasta进行数据预处理
  • 2: 使用DIAMOND进行all-vs-all blast
  • 3: 使用MySQL数据库寻找配对
  • 4: 使用MCL处理配对

分析流程

我们以一个具体的案例说明下。新建一个文件夹，例如说test， 里面是收集好的物种氨基酸序列。

1
2
$ ls
Alyrata.faa  Athaliana.faa  Chirsuta.faa

由于Docker版本的OrthoMCL使用DIAMOND进行blastp比对，因此一定要保证你的氨基酸序列中没有”.”，否则会报错。可以用seqkit grep过滤不合格的序列。

1
seqkit grep -s -vrp '"\."' input.fa > output.faa

之后创建一个fasta_table，放在test目录下。该文件分为两列，第一列是文件名，第二列是缩写。

1
2
3
Alyrata.faa	Aly
Athaliana.faa	Ath
Chirsuta.faa	Chi

准备配置文件orthomcl.config，同样放在test目录下

1
2
3
4
5
6
7
8
9
10
11
12
dbVendor=mysql
dbConnectString=dbi:mysql:orthomcl:mysql_local_infile=1
dbLogin=root
dbPassword=PAssw0rd
similarSequencesTable=SimilarSequences
orthologTable=Ortholog
inParalogTable=InParalog
coOrthologTable=CoOrtholog
interTaxonMatchView=InterTaxonMatch
percentMatchCutoff=50
evalueExponentCutoff=-5
oracleIndexTblSpc=NONE

在test同级目录下运行如下命令

1
docker run --privileged=true  --rm -it --volume $PWD/test:/outdir:rw jasonkwan/orthomcl_docker:latest bash

就会进入Docker的交互命令行，运行run_orthomcl.py

1
2
cd outdir
run_orthomcl.py --table_path fasta_table --config_file orthomcl.config --processors 32 &

假如不希望进入交互命令行，那么需要按照下面的方法进行运行

1
docker run --privileged=true  --rm --volume $PWD/test:/outdir:rw jasonkwan/orthomcl_docker:latest /bin/bash -c "/tmp/.runconfig.sh && run_orthomcl.py --table_path /outdir/fasta_table --config_file /outdir/orthomcl.config --processors 32 "

最终输出结果是groups.txt。下一个问题是，当你有了groups.txt后，下面能做什么分析呢？

按照对基因组学文章的整理，基本上就是下面两个

  • 使用单拷贝基因构建系统发育树。
  • 使用CAFE进行基因家族扩张收缩分析

我正在找一篇文章尝试重现这个流程。

参考资料

生信软件工具箱

这可能是最全最好的BLAST教程

Basic local alignment search tool (BLAST)

包括：blastn, blastp, blastx, tblastn, tblastx等. 使用conda安装即可。

1
2
3
conda install -c bioconda blast
# blast安装perl模块的方法
conda install perl-digest-md5

BLAST的主要理念

  • Search may take place in nucleotide and/or protein space or translated spaces where nucleotides are translated into proteins.
  • Searches may implement search “strategies”: optimizations to a certain task. Different search strategies will return different alignments.
  • Searches use alignments that rely on scoring matrices
  • Searches may be customized with many additional parameters. BLAST has many subtle functions that most users never need.

本地BLAST的基本步骤

  1. 用makeblastdb为BLAST提供数据库
  2. 选择blast工具，blastn,blastp
  3. 运行工具，有必要的还可以对输出结果进行修饰

第一步：建立检索所需数据库

BLAST数据库分为两类，核酸数据库和氨基酸数据库，可以用makeblastbd创建。可以用help参数简单看下说明。

1
2
3
4
5
6
7
8
9
10
$ makeblastdb -help
USAGE
  makeblastdb [-h] [-help] [-in input_file] [-input_type type]
    -dbtype molecule_type [-title database_title] [-parse_seqids]
    [-hash_index] [-mask_data mask_data_files] [-mask_id mask_algo_ids]
    [-mask_desc mask_algo_descriptions] [-gi_mask]
    [-gi_mask_name gi_based_mask_names] [-out database_name]
    [-max_file_sz number_of_bytes] [-logfile File_Name] [-taxid TaxID]
    [-taxid_map TaxIDMapFile] [-version]
-dbtype <String, `nucl', `prot'>

具体以拟南芥基因组作为案例，介绍使用方法：
注： 拟南芥的基因组可以在TAIR上下在，也可在ensemblgenomes下载。后者还可以下载其他植物的基因组

1
2
3
4
5
6
7
8
9
10
11
# 下载基因组
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-36/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
gzip -d Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
# 构建核酸BLAST数据库
makeblastdb -in Arabidopsis_thaliana.TAIR10.dna.toplevel.fa -dbtype nucl -out TAIR10 -parse_seqids

# 下载拟南芥protein数据
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-36/fasta/arabidopsis_thaliana/pep/Arabidopsis_thaliana.TAIR10.pep.all.fa.gz
# 构建蛋白BLAST数据库
gzip -dArabidopsis_thaliana.TAIR10.pep.all.fa.gz
makeblastdb -in  Arabidopsis_thaliana.TAIR10.pep.all.fa -dbtype prot -out TAIR10 -parse_seqids

如果你从NCBI或者其他渠道下载了格式化过的数据库，那么可以用blastdbcmd去检索blast数据库，参数很多，常用就如下几个

  • db string : string表示数据库所在路径
  • dbtype string,: string在（guess, nucl, prot)中选择一个
  • 检索相关参数
    • -entry all 或 555, AC147927 或 gnl|dbname|tag’
    • -entry_batch 提供一个包含多个检索关键字的文件
    • -info 数据库基本信息
  • 输出格式 -outfmt %f %s %a %g …默认是%f
  • out 输出文件
  • show_blastdb_search_path: blast检索数据库路径

使用案例

1
2
3
4
5
6
# 查看信息
blastdbcmd -db TAIR10 -dbtype nucl -info
# 所有数据
blastdbcmd -db TAIR10 -dbtype nucl -entry all | head
# 具体关键字，如GI号
blastdbcmd -db TAIR10 -dbtype nucl -entry 3 | head

还有其他有意思的参数，可以看帮助文件了解

可选：BLAST安装和更新nr和nt库

安装nt/nr库需要先进行环境变量配置，在家目录下新建一个.ncbirc配置文件，然后添加如下内容

1
2
3
4
5
6
7
8
9
10
11
12
13
; 开始配置BLAST
[BLAST]
; 声明BLAST数据库安装位置
BLASTDB=/home/xzg/Database/blast
; Specifies the data sources to use for automatic resolution
; for sequence identifiers
DATA_LOADERS=blastdb
; 蛋白序列数据库存本地位置
BLASTDB_PROT_DATA_LOADER=/home/xzg/Database/blast/nr
; 核酸数据库本地存放位置
BLASTDB_NUCL_DATA_LOADER=/home/xzg/Database/blast/nt
[WINDOW_MASKER]
WINDOW_MASKER_PATH=/home/xzg/Database/blast/windowmasker

配置好之后，使用BLAST+自带的update_blastdb.pl脚本下载nr和nt等库文件（不建议下载序列文件，一是因为后者文件更大，二是因为可以从库文件中提取序列blastdbcmd -db nr -dbtype prot -entry all -outfmt “%f” -out nr.fa ，最主要是建库需要花费很长时间），直接运行下列命令即可自动下载。

1
nohup time update_blastdb.pl nt nr > log &

如果你不像通过update_blastdb.pl下载nr和nt等库文件，也可以是从ncbi上直接下载一系列nt/nr.xx.tar.gz文件，然后解压缩即可，后续还可以用update_blastdb.pl进行数据更新。

报错： 使用update_blastdb.pl更新和下载数据库时候出现模块未安装的问题。解决方法，首先用conda安装对应的模块，然后修改update_blastdb.pl的第一行，即shebang部分，以conda的perl替换，或者按照如下方法执行。

1
perl `which update_blastdb.pl`

下载过程中请确保网络状态良好，否则会出现Downloading nt.00.tar.gz...Unable to close datastream报错。

第二步：选择blast工具

根据不同的需求，比如说你用的序列是氨基酸还是核苷酸，你要查找的数据是核甘酸还是氨基酸，选择合适的blast工具。不同需求的对应关系可以见下图（来自biostars handbook）

BLAST工具

不同工具的应用范围虽然不同，但是基本参数都是一致的，学会blastn，基本上其他诸如blastp,blastx也都会了。

blastn的使用参数很多 blastn [-h] ，但是比较常用是如下几个

  • -db : 数据库在本地的位置，或者是NCBI上数据库的类型,如 -db nr

BLAST库

  • -query: 检索文件

  • -query_loc : 指定检索的位置

  • -strand： 搜索正义链还是反义链，还是都要

  • out : 输出文件

  • -remote： 可以用NCBI的远程数据库， 一般与 -db nr

  • -evalue 科学计数法，比如说1e3，定义期望值阈值。E值表明在随机的情况下，其它序列与目标序列相似度要大于这条显示的序列的可能性。 与S值有关，S值表示两序列的同源性，分值越高表明它们之间相似的程度越大

    E值总结： 1.E值适合于有一定长度，而且复杂度不能太低的序列。2. 当E值小于10-5时，表明两序列有较高的同源性，而不是因为计算错误。3. 当E值小于10-6时，表时两序列的同源性非常高，几乎没有必要再做确认。

  • 与联配计分相关参数： -gapopen，打开gap的代价；-gapextend， gap延伸的代价；-penalty：核酸错配的惩罚； -reward， 核酸正确匹配的奖励；

  • 结果过滤：-perc_identity, 根据相似度

注 BLAST还提供一个task参数，感觉很有用的样子，好像会针对任务进行优化速度。

第三步：运行blast，调整输出格式。

我随机找了一段序列进行检索

1
2
echo '>test' > query.fa
echo TGAAAGCAAGAAGAGCGTTTGGTGGTTTCTTAACAAATCATTGCAACTCCACAAGGCGCCTGTAATAGACAGCTTGTGCATGGAACTTGGTCCACAGTGCCCTACCACTGATGATGTTGATATCGGAAAGTGGGTTGCAAAAGCTGTTGATTGTTTGGTGATGACGCTAACAATCAAGCTCCTCTGGT >> query.fa

用的是blastn 检索核酸数据库。最简单的用法就是提供数据库所在位置和需要检索的序列文件

1
2
3
4
5
6
blastn -db BLAST/TAIR10 -query query.fa
# 还可以指定检索序列的位置
blastn -db BLAST/TAIR10 -query query.fa  -query_loc 20-100
# 或者使用远程数据库
blastn -db nr -remote -query query.fa
blastn -db nt -remote -query query.fa

以上是默认输出，blast的-outfmt选项提供个性化的选择。一共有18个选择，默认是0。
0 = Pairwise,
1 = Query-anchored showing identities,
2 = Query-anchored no identities,
3 = Flat query-anchored showing identities,
4 = Flat query-anchored no identities,
5 = BLAST XML,
6 = Tabular,
7 = Tabular with comment lines,
8 = Seqalign (Text ASN.1),
9 = Seqalign (Binary ASN.1),
10 = Comma-separated values,
11 = BLAST archive (ASN.1),
12 = Seqalign (JSON),
13 = Multiple-file BLAST JSON,
14 = Multiple-file BLAST XML2,
15 = Single-file BLAST JSON,
16 = Single-file BLAST XML2,
17 = Sequence Alignment/Map (SAM),
18 = Organism Report
其中6，7，10，17可以自定输出格式。默认是

qaccver saccver pident length mismatch gapopen qstart qend sstart send evalue bitscore

简写 含义
qaccver 查询的AC版本（与此类似的还有qseqid,qgi,qacc,与序列命名有关）
saccver 目标的AC版本（于此类似的还有sseqid,sallseqid,sgi,sacc,sallacc,也是序列命名相关)
pident 完全匹配百分比 （响应的nident则是匹配数）
length 联配长度(另外slen表示查询序列总长度，qlen表示目标序列总长度)
mismatch 错配数目
gapopen gap的数目
qstart 查询序列起始
qstart 查询序列结束
sstart 目标序列起始
send 目标序列结束
evalue 期望值
bitscore Bit得分
score 原始得分
AC： accession

以格式7为实例进行输出，并且对在线数据库进行检索 

1
2
blastn -task blastn -remote -db nr -query query.fa  -outfmt 7 -out query.txt
head -n 15 query.txt

BLAST BEST HIT

如果想输入序列，增加对应的格式qseq, sseq

1
blastn -query query.fa -db nr -outfmt ' 6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send qseq sseq evalue bitscore'

对已有序列进行注释时常见的best hit only模式命令行

1
2
3
4
5
6
7
blastn -query gene.fa -out gene.blast.txt -task megablast -db nt -num_threads 12 -evalue 1e-10 -best_hit_score_edge 0.05 -best_hit_overhang 0.25 -outfmt "7 std stitle" -perc_identity 50 -max_target_seqs 1
# 参数详解
-task megablast : 任务执行模式，可选有'blastn' 'blastn-short' 'dc-megablast' 'megablast' 'rmblastn'
-best_hit_score_edge 0.05 :  Best Hit 算法的边界值，取值范围为0到0.5，系统推荐0.1
-best_hit_overhang 0.25 : Best Hit 算法的阈值，取值范围为0到0.5，系统推荐0.1
-perc_identity 50 : 相似度大于50
-max_target_seqs 1 : 最多保留多少个联配

仅仅看参数依旧无用，还需要知道BLAST的Best-Hits的过滤算法。假设一个序列存在两个match结果,A和B，无论A还是B，他们的HSP（High-scoring Segment Pair, 没有gap时的最高联配得分）一定要高于best_hit_overhang，否则被过滤。如果满足下列条件则保留A

  • evalue(A) >= evalue (B)
  • score(A)/length(A) < (1.0-score_edge)*score(B)/length(B)

搭建网页BLAST

曾经的BLAST安装后提供wwwblast用于构建本地的BLAST网页工具，但是BLAST+没有提供这个工具，好在BLAST足够出名，也就有人给它开发网页版工具。如viroBLAST和Sequenceserver, 目前来看似乎后者更受人欢迎。

有root安装起来非常的容易

1
2
sudo apt install ruby ncbi-blast+ ruby-dev rubygems-integration npm
sudo gem install sequenceserver

数据库准备,前面步骤已经下载了拟南芥基因组的FASTA格式数据

1
sequenceserver -d /到/之前/建立/BLAST/文件路径

然后就可以打开浏览器输入IP:端口号使用了。

Sequenceserver高级用法

开机启动：

新建一个/etc/systemd/system/sequenceserver.service文件，添加如下内容。注意修改ExecStart.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
[Unit]
Description=SequenceServer server daemon
Documentation="file://sequenceserver --help" "http://sequenceserver.com/doc"
After=network.target

[Service]
Type=simple
User=seqservuser
ExecStart=/path/to/bin/sequenceserver -c /path/to/sequenceserver.conf
KillMode=process
Restart=on-failure
RestartSec=42s
RestartPreventExitStatus=255

[Install]
WantedBy=multi-user.target

然后重新加载systemctl

1
2
3
4
5
6
## let systemd know about changed files
sudo systemctl daemon-reload
## enable service for automatic start on boot
systemctl enable sequenceserver.service
## start service immediately
systemctl start sequenceserver.service

nginx反向代理：我承认没有基本的nginx的知识根本搞不定这一步，所以我建议组内使用就不要折腾这个。简单的说就是在nginx的配置文件的server部分添加如下内容即可。

1
2
3
4
5
6
7
location /sequenceserver/ {
    root /home/priyam/sequenceserver/public/dist;
    proxy_pass http://localhost:4567/;
    proxy_intercept_errors on;
    proxy_connect_timeout 8;
    proxy_read_timeout 180;
}

假如你有docker的话

1
2
3
4
# ubuntu
docker run --rm -itp 4567:4567 -v /path-to-database-dir:/db wurmlab/sequenceserver:1.0.11
# centos
docker run --privileged --rm -itp 4567:4567 -v /path-to-database-dir:/db wurmlab/sequenceserver:1.0.11

注意，你得在宿主机器上 /path-to-database-dir建好blast索引，或者用docker exec -it docker容器名 /bin/bash 进入到容器里，用sequenceserver -d db &运行

参考资料：http://www.sequenceserver.com/doc/

术语列表

引自BLAST Glossary

参考资料

生信软件工具箱

「JCVI教程」使用JCVI根据CDS的BLAST结果绘制共线性图

这是唐海宝老师GitHub上的JCVI工具的非官方说明书。
该工具集的功能非常多，但是教程资料目前看起来并不多，因此为了能让更多人用上那么好用的工具，我就一边探索，一边写教程

这一篇文章教大家如何利用JCVI里面的工具绘制点图，展现两个物种之间的共线性关系。

在分析之前，你需要从PhytozomeV11 下载A.thaliana和Alyrata的CDS序列，保证文件夹里有如下内容

1
2
Alyrata_384_v2.1.cds.fa.gz        Athaliana_167_TAIR10.cds.fa.gz
Alyrata_384_v2.1.gene.gff3.gz  Athaliana_167_TAIR10.gene.gff3.gz

准备最长CDS和BED文件

我们在做CDS相互比对的时候只需要有每个基因最长的转录本即可，有两种方法可以实现

方法1：自己写脚本

我用我写的一个脚本get_the_longest_transcripts.py提取每个基因的最长转录本,见 基因组共线性工具MCScanX使用说明

1
2
zcat Alyrata_384_v2.1.gene.gff3.gz | python ~/scripts/python/get_the_longest_transcripts.py > aly_lst_gene.txt
zcat Athaliana_167_TAIR10.gene.gff3.gz | python ~/scripts/python/get_the_longest_transcripts.py  > ath_lst_gene.txt

其中xxx_lst_gene.txt的格式如下, 第一列是基因名，第二列是mRNA编号，后面几列是位置信息。

1
2
3
4
5
6
7
8
9
10
11
$ head ath_lst_gene.txt
AT4G19470.TAIR10	AT4G19470.1.TAIR10	Chr4	10612993	10614339	-
AT5G43860.TAIR10	AT5G43860.1.TAIR10	Chr5	17630450	17632312	+
AT1G68650.TAIR10	AT1G68650.1.TAIR10	Chr1	25775741	25777874	+
AT1G28050.TAIR10	AT1G28050.1.TAIR10	Chr1	9775528	9777810	-
AT3G59880.TAIR10	AT3G59880.1.TAIR10	Chr3	22120969	22121700	+
AT1G22030.TAIR10	AT1G22030.1.TAIR10	Chr1	7759164	7760556	-
AT5G24330.TAIR10	AT5G24330.1.TAIR10	Chr5	8295147	8297068	-
AT5G43990.TAIR10	AT5G43990.2.TAIR10	Chr5	17697889	17702005	+
AT1G11410.TAIR10	AT1G11410.1.TAIR10	Chr1	3841286	3844432	+
AT4G32890.TAIR10	AT4G32890.1.TAIR10	Chr4	15875470	15876762	+

由于基因名和mRNA编号里有在提取CDS不需要的内容，因此要进行删除

1
2
sed -i 's/\.v2\.1//g' aly_lst_gene.txt
sed -i 's/\.TAIR10//g' ath_lst_gene.txt

之后我们就可以根据第二列进行提取CDS

1
2
seqkit grep -f  <(cut -f 2 ath_lst_gene.txt ) Athaliana_167_TAIR10.cds.fa.gz > ath.cds
seqkit grep -f  <(cut -f 2 aly_lst_gene.txt ) Alyrata_384_v2.1.cds.fa.gz > aly.cds

提取的CDS编号里面也有一些不需要的内容，所以也要删除

1
2
sed -i 's/\.t.*//' aly.cds
sed -i 's/\..*//' ath.cds

此外还需要基因的位置信息的bed文件

1
2
awk '{print $3"\t"$4"\t"$5"\t"$1"\t0\t"$6}' ath_lst_gene.txt | sort -k4,4V > ath.bed
awk '{print $3"\t"$4"\t"$5"\t"$1"\t0\t"$6}' aly_lst_gene.txt | sort -k4,4V > aly.bed

基于JCVI工具集

当然也可以参考「JCVI教程」如何基于编码序列或蛋白序列进行共线性分析来提取bed和cds序列，不需要用到我写的脚本。

1
2
python -m jcvi.formats.gff bed --type=mRNA --key=Name Athaliana_167_TAIR10.gene.gff3.gz  > ath.bed
python -m jcvi.formats.gff bed --type=mRNA --key=Name Alyrata_384_v2.1.gene.gff3.gz > aly.bed

对bed文件中的基因进行去重

1
2
python -m jcvi.formats.bed uniq ath.bed
python -m jcvi.formats.bed uniq aly.bed

这一步会得到aly.uniq.bed和ath.uniq.bed, 我们将其覆盖原文件

1
2
mv ath.uniq.bed ath.bed
mv aly.uniq.bed aly.bed

根据bed文件提取cds里的序列

1
2
3
4
# Athaliana
seqkit grep -f <(cut -f 4 ath.bed ) Athaliana_167_TAIR10.cds.fa.gz | seqkit seq -i > ath.cds
# Alyrata
seqkit grep -f <(cut -f 4 aly.bed )  Alyrata_384_v2.1.cds.fa.gz | seqkit seq -i  > aly.cds

无论是哪种方法，请保证最后有以下四个文件

1
2
$ ls ???.???
aly.bed  aly.cds  ath.bed  ath.cds

BLAST比对

相对于上一步，这一步其实非常简单了

1
2
makeblastdb -in ath.cds -out db/ath -dbtype nucl
blastn -num_threads 20  -query aly.cds -db db/ath -outfmt 6 -evalue 1e-5 -num_alignments 5 > aly_ath.blast

用jcvi.compara.blastfilter 对结果进行过滤

1
python -m jcvi.compara.blastfilter --no_strip_names aly_ath.blast --sbed ath.bed --qbed aly.bed

运行过程中有如下输出信息

1
2
3
4
5
6
7
8
19:59:15 [base] Load file `aly.bed`
19:59:16 [base] Load file `ath.bed`
19:59:16 [blastfilter] Load BLAST file `aly_ath.blast` (total 49887 lines)
19:59:16 [base] Load file `aly_ath.blast`
19:59:16 [blastfilter] running the cscore filter (cscore>=0.70) ..
19:59:16 [blastfilter] after filter (42023->26531) ..
19:59:16 [blastfilter] running the local dups filter (tandem_Nmax=10) ..
19:59:16 [blastfilter] after filter (26531->24242) ..

最后输出aly_ath.blast.filtered用于做图

1
python -m jcvi.graphics.blastplot aly_ath.blast.filtered --sbed ath.bed --qbed aly.bed

最后点图如下

点图

生信软件工具箱

「JCVI教程」使用JCVI进行基因组共线性分析(中)

在「JCVI教程」使用JCVI进行基因组共线性分析(上)还是有一个尴尬的事情，就是只用到两个物种，不能展示出JCVI画图的方便之处，因此这里参考https://github.com/tanghaibao/jcvi/wiki/MCscan-(Python-version)的分析，只不过画图部分拓展下思路。

首先要运行如下代码获取目的数据

1
2
3
4
5
6
7
8
9
10
11
12
13
python -m jcvi.apps.fetch phytozome Vvinifera,Ppersica,Tcacao
python -m jcvi.formats.gff bed --type=mRNA --key=Name Vvinifera_145_gene.gff3.gz -o grape.bed
python -m jcvi.formats.gff bed --type=mRNA --key=Name Ppersica_139_gene.gff3.gz -o peach.bed
python -m jcvi.formats.gff bed --type=mRNA --key=Name Tcacao_233_gene.gff3.gz -o cacao.bed
python -m jcvi.formats.fasta format --sep="|" Vvinifera_145_cds.fa.gz grape.cds
python -m jcvi.formats.fasta format --sep="|" Ppersica_139_cds.fa.gz peach.cds
python -m jcvi.formats.fasta format --sep="|" Tcacao_233_cds.fa.gz cacao.cds
# find ortholog
python -m jcvi.compara.catalog ortholog grape peach --cscore=.99
python -m jcvi.compara.catalog ortholog peach cacao --cscore=.99
# build .simpe
python -m jcvi.compara.synteny screen --minspan=30 --simple grape.peach.anchors grape.peach.anchors.new 
python -m jcvi.compara.synteny screen --minspan=30 --simple peach.cacao.anchors peach.cacao.anchors.new

然后按照教程的配置文件进行画图

layout文件内容如下

1
2
3
4
5
6
7
# y, xstart, xend, rotation, color, label, va,  bed
 .7,     .2,    .4,      45,      , Grape, top, grape.bed
 .5,     .2,    .8,       0,      , Peach, top, peach.bed
 .7,     .4,    .8,     -45,      , Cacao, bottom, cacao.bed
# edges
e, 0, 1, grape.peach.anchors.simple
e, 1, 2, peach.cacao.anchors.simple

seqids文件内容如下

1
2
3
chr1,chr2,chr3,chr4,chr5,chr6,chr7,chr8,chr9,chr10,chr11,chr12,chr13,chr14,chr15,chr16,chr17,chr18,chr19
scaffold_1,scaffold_2,scaffold_3,scaffold_4,scaffold_5,scaffold_6,scaffold_7,scaffold_8
scaffold_1,scaffold_2,scaffold_3,scaffold_4,scaffold_5,scaffold_6,scaffold_7,scaffold_8,scaffold_9,scaffold_10r

运行python -m jcvi.graphics.karyotype seqids layout会得到如下结果

三个物种的共线性图

我就在思考一个问题如何让他形成一个三角形。经过一波三角运算和不断尝试，我定义了如下的layout

1
2
3
4
5
6
7
# y, xstart, xend, rotation, color, label, va,  bed
 .5,      0.025,    0.625,      60,      , Grape, top, grape.bed
 .2,      0.2,    .8,       0,      , Peach, top, peach.bed
 .5,     0.375,    0.975,     -60,      , Cacao, top, cacao.bed
# edges
e, 0, 1, grape.peach.anchors.simple
e, 1, 2, peach.cacao.anchors.simple

那么效果怎么样呢？运行python -m jcvi.graphics.karyotype seqids layout吧

三角形

当然这里只展示了，grape和peach, peach和cacao之间的共线性，我又想着能不能加上grape和caocao呢？我尝试着运行下面的代码，

1
2
python -m jcvi.compara.catalog ortholog grape cacao --cscore=.99
python -m jcvi.compara.synteny screen --minspan=30 --simple grape.cacao.anchors grape.cacao.anchors.new

并继续修改了layout

1
2
3
4
5
6
7
8
# y, xstart, xend, rotation, color, label, va,  bed
 .5,      0.025,    0.625,      60,      , Grape, top, grape.bed
 .2,      0.2,    .8,       0,      , Peach, top, peach.bed
 .5,     0.375,    0.975,     -60,      , Cacao, top, cacao.bed
# edges
e, 0, 1, grape.peach.anchors.simple
e, 1, 2, peach.cacao.anchors.simple
e, 0, 2, grape.cacao.anchors.simple

运行python -m jcvi.graphics.karyotype seqids layout会得到如下结果

awesome

我觉得给我一点时间，我也能用JCVI画出下面的图了

amazing

图片来自于https://science.sciencemag.org/content/345/6199/950

生信软件工具箱

「JCVI教程」使用JCVI进行基因组共线性分析(上)

本教程借鉴https://github.com/tanghaibao/jcvi/wiki/MCscan-(Python-version).

我们先从http://plants.ensembl.org/index.html选择两个物种做分析, 这里选择的就是前两个物种，也就是拟南芥和水稻（得亏没有小麦和玉米）

选择物种

我们下载它的GFF文件，cdna序列和蛋白序列

1
2
3
4
5
6
7
8
#Athaliana
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/fasta/arabidopsis_thaliana/cdna/Arabidopsis_thaliana.TAIR10.cdna.all.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/fasta/arabidopsis_thaliana/pep/Arabidopsis_thaliana.TAIR10.pep.all.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.44.gff3.gz
#Osativa
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/fasta/oryza_sativa/cdna/Oryza_sativa.IRGSP-1.0.cdna.all.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/fasta/oryza_sativa/pep/Oryza_sativa.IRGSP-1.0.pep.all.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/gff3/oryza_sativa/Oryza_sativa.IRGSP-1.0.44.gff3.gz

保证要有6个文件以便下游分析

1
2
3
$ ls
Arabidopsis_thaliana.TAIR10.44.gff3.gz      Arabidopsis_thaliana.TAIR10.pep.all.fa.gz  Oryza_sativa.IRGSP-1.0.cdna.all.fa.gz
Arabidopsis_thaliana.TAIR10.cdna.all.fa.gz  Oryza_sativa.IRGSP-1.0.44.gff3.gz          Oryza_sativa.IRGSP-1.0.pep.all.fa.gz

我们分析只需要用到每个基因最长的转录本就行，之前我用的是自己写的脚本，但其实我发现jcvi其实可以做到这件事情

先将gff转成bed格式，

1
2
python -m jcvi.formats.gff bed --type=mRNA --key=transcript_id Arabidopsis_thaliana.TAIR10.44.gff3.gz > ath.bed
python -m jcvi.formats.gff bed --type=mRNA --key=transcript_id Oryza_sativa.IRGSP-1.0.44.gff3.gz > osa.bed

然后将bed进行去重复

1
2
python -m jcvi.formats.bed uniq ath.bed
python -m jcvi.formats.bed uniq osa.bed

最后我们得到了ath.uniq.bed和osa.uniq.bed, 根据bed文件第4列就可以用于提取cds序列和蛋白序列。

1
2
3
4
5
6
# Athaliana
seqkit grep -f <(cut -f 4 ath.uniq.bed ) Arabidopsis_thaliana.TAIR10.cdna.all.fa.gz | seqkit seq -i > ath.cds
seqkit grep -f <(cut -f 4 ath.uniq.bed ) Arabidopsis_thaliana.TAIR10.pep.all.fa.gz | seqkit seq -i > ath.pep 
# Osativa
seqkit grep -f <(cut -f 4 osa.uniq.bed )  Oryza_sativa.IRGSP-1.0.cdna.all.fa.gz | seqkit seq -i  > osa.cds
seqkit grep -f <(cut -f 4 osa.uniq.bed ) Oryza_sativa.IRGSP-1.0.pep.all.fa.gz | seqkit seq -i  > osa.pep

这里用到的seqkit建议学习，非常好用

下面使用python -m jcvi.compara.catalog ortholog进行共线性分析，这是一个非常行云流水的过程（除非你报错）

新建一个文件夹，方便在报错的时候，把全部都给删了，

1
2
3
4
5
mkdir -p cds && cd cds
ln -s ../ath.cds ath.cds
ln -s ../ath.uniq.bed ath.bed
ln -s ../osa.cds osa.cds
ln -s ../osa.uniq.bed osa.bed

运行代码

1
python -m jcvi.compara.catalog ortholog --no_strip_names ath osa

输出结果如下

1
2
$ ls ath.osa.*
ath.osa.anchors  ath.osa.last  ath.osa.last.filtered  ath.osa.lifted.anchors  ath.osa.pdf

其中我们最感兴趣都是pdf结果，不出意外没啥共线性。

共线性结果

我们还可以用蛋白序列做共线性分析

1
2
3
4
5
6
# 在之前输出cds,pep都文件夹操作
mkdir -p pep && cd pep
ln -s ../ath.pep ath.pep
ln -s ../ath.uniq.bed ath.bed
ln -s ../osa.pep osa.pep
ln -s ../osa.uniq.bed osa.bed

运行代码

1
python -m jcvi.compara.catalog ortholog --dbtype prot --no_strip_names ath osa

我之前以为他不可以基于蛋白序列分析，幸亏有人提醒。

蛋白共线性

你会发现这是一个自动化分析流程，我们只是提供了4个文件，它就完成了一些列事情。它生成的文件里除了PDF外，其他还有啥用呢？

  • ath.osa.last: 基于LAST的比对结果
  • ath.osa.last.filtered: LAST的比对结果过滤串联重复和低分比对
  • ath.osa.anchors: 高质量的共线性区块
  • ath.osa.lifted.anchors:增加了额外的锚点，形成最终的共线性区块

anchors文件特别有用，之后会写一篇介绍如何利用他进行可视化，这里介绍它的格式。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
###
AT1G28395.5     Os01t0238800-02 66
AT1G28440.1     Os01t0239700-02 1360
AT1G28480.1     Os01t0241400-01 136
AT1G28510.1     Os01t0242300-01 241
###
AT1G11100.3     Os01t0779400-01 943
AT1G11125.1     Os01t0779800-01 52
AT1G11160.2     Os01t0780400-02 535
AT1G11180.1     Os01t0780500-01 483
AT1G11330.2     Os01t0784700-00 742
AT1G11360.1     Os01t0783500-01 305
AT1G11540.2     Os01t0786800-01 422
AT1G11570.3     Os01t0788200-01 162
AT1G11580.2     Os01t0788400-01 550
AT1G11630.1     Os01t0793200-01 321

每个共线性区块以###进行分隔, 第一列是检索的基因，第二列是被检索的基因，第三列则是两个序列的BLAST的bit score,值越大可靠性越高。

用水稻和拟南芥进行了比较之后，发现后面基本上也没啥可以分析了。因此下面基于「JCVI教程」如何基于物种的CDS的blast结果绘制点图(dotplot)得到的cds和bed文件进行分析。

之前已经得到了如下四个文件

1
2
ls ???.???
aly.bed  aly.cds  ath.bed  ath.cds

所以我们只要运行

1
python -m jcvi.compara.catalog ortholog --no_strip_names aly ath

就得到了一个非常好看的点图

共线性点图

我们可以发现，都作为Arabidopsis属的两个物种，他们之间存在很高的同源性，并且同源区比例是1:1,

共线区域比例

这其实和2011年的Nature Genetics上Alyrata的文章的结果是相似的，只不过他不是用点图进行展示

Nature Genetics

我们也可以用JCVI的画图模块实现这种效果，只不过还需要一点额外操作，创建如下三个文件

  • seqids: 需要展现哪些序列
  • layout: 不同物种的在图上的位置
  • .simple: 从.anchors文件创建的更简化格式

第一步，创建.simple文件

1
python -m jcvi.compara.synteny screen --minspan=30 --simple aly.ath.anchors aly.ath.anchors.new

第二步, 创建seqid文件，非常简单，就是需要展示的scaffold或染色体的编号

1
2
scaffold_1,scaffold_2,scaffold_3,scaffold_4,scaffold_5,scaffold_6,scaffold_7,scaffold_8
Chr1,Chr2,Chr3,Chr4,Chr5

第二步，创建layout文件，用于设置绘制的一些选项。

1
2
3
4
5
# y, xstart, xend, rotation, color, label, va,  bed
 .6,     .2,    .8,       0,      , Alyrata, top, aly.bed
 .4,     .2,    .8,       0,      , Athaliana, top, ath.bed
# edges
e, 0, 1, aly.ath.anchors.simple

注意, #edges下的每一行开头都不能有空格

最后运行下面的命令，会得到一个karyotype.pdf

1
python -m jcvi.graphics.karyotype seqids layout

染色体共线性图

如何让这个图垂直呢？（导入AI里就好了）