基因家族鉴定流程（一）数据准备

NCBI高级搜索以及序列下载

参照https://www.ncbi.nlm.nih.gov/books/NBK49540/以及https://www.ncbi.nlm.nih.gov/books/NBK3837/所展示的高级搜索方法。
以人类Globin为例，如果直接以"human"[Organism] AND globin为关键词进行搜索，将会得到如下结果

fig1.png

可以发现一共有926个结果，其中包括了很多同一基因的不同转录本，以及来自不同数据库的相同基因对应的转录本，以及并非是globin的蛋白，这是因为如果你仅仅以globin为关键词，那么会搜索text、Journal、title等多个区段，有些蛋白可能是和globin有相关的作用，而记录在这些区段当中，从而被检索到。我们可以设置关键词，只搜索refseq结果，因为refseq结果的质量是要比genbank这些数据库的结果要好很多的。

fig2.png

可以看到，现在结果数量缩小了很多，但仍然存在不是globin的蛋白以及冗余转录本。
如果你搜索的物种是人类、小鼠这样的物种，那么恭喜你，你可以用Refseq Select[filter]来得到refseq代表转录本。如下图所示：

fig3.png

目前Refseq Select只支持小鼠、人类以及原核生物，NCBI官方未来计划包含更多的模式物种以及额外的真核生物物种，详见：
https://www.ncbi.nlm.nih.gov/refseq/refseq_select
如果你想检索其他物种，将最长转录本作为代表，那就有点麻烦了，因为光从protein这个搜索框是没办法挑选出最长或者以其他方式作为代表的转录本所对应的蛋白的，命名的不同以及各种原因导致很多时候我们需要手动筛选，笨办法就是通过gene搜索框搜索相关基因，然后可以看到基因相关的不同转录本，选择最长的一个个下载，这样会比较准确但是比较慢，有的时候我还是宁愿花很多时间去慢慢手工筛选的，确保不漏掉任何一个蛋白，不多任何一个其他不相干的蛋白。这里我们先继续以人类为例。
Tips：如果你明确知道Gene_symbolID名，就可以直接用AND xxx[GENE]搜索，会得到有冗余转录本但都是相关基因的蛋白，但是由于历史原因（命名的多种多样）等等，这种方法检索虽然准确但是效率极低，如果你确实只需要找某个亚家族子成员，这种方法会比较合适。
目前我们希望得到的检索结果的名称，是仅仅包含globin这一关键词的条目（当然有些蛋白可能它是这个家族的成员，但是DEFINITION字段的命名并非和其他成员一样，因此也要考虑很多方面，需要你对这个家族足够了解），我们可以先下载summry

fig4.png

然后进入linux环境，处理一下。

# 忽略大小写，搜索关键词，globin，然后输出RefSeq ID，注意有的虽然包含了globin，但是仍然可能不是这个家族的成员，因此还是最好人工审核一下
gawk 'BEGIN{IGNORECASE=1;RS="[0-9][.] ";FS="\n"}$1~/globin/{print $3}' protein_result.txt | cut -d " " -f 1 > id_list.txt

这样我们得到了带有globin检索词的相关的蛋白序列ID号，

fig5.png

接下来我们可以用NCBI官方的下载工具Edirect下载序列
注：shell用的是zsh而不是bash

# 如果没有可以用conda安装
conda install -c bioconda entrez-direct
# 下载序列
for id (${(f)"$(<id_list.txt)"}) { 
    esearch -db protein -query $id | 
    efetch -format fasta >> globin.faa
}

基因组&&目标cds获取

一般常用的基因组下载来源就是ensemble和NCBI genome，当然还有很多其他的数据库。
这里以小鼠的基因组为例，在genome搜索框搜索点进相应条目即可，这里介绍另一种下载方式，我们右键单击GFF，点击复制链接地址

fig6.png

然后我们粘贴到地址栏，去掉选中的这一部分，即可打开新大陆，即ftp下载地址

fig7.png

fig8.png

想要知道各自对应的是什么信息，可以点击README.txt查看，我们可以点击gff.gz和fna.gz下载下来相应的基因组和gff注释然后提取CDS（用gffread或者写个简单的脚本都可以），或者直接
下载translated.cds.faa.gz，但是这里如果你直接下的cds，那么会包含冗余的转录本序列，所以最好的方法还是下载基因组和gff，然后写个脚本提取最长转录本，这个脚本写法也很简单网上也有很多人的perl/python脚本，可以参考https://www.jianshu.com/p/f173adf932ae这位大哥的脚本。这里我为了简便起见直接下cds。

# 获取翻译后的cds
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/635/GCF_000001635.27_GRCm39/GCF_000001635.27_GRCm39_translated_cds.faa.gz
# 解压,得到相应序列
gunzip GCF_000001635.27_GRCm39_translated_cds.faa.gz

Uniprot介绍及序列检索下载

介绍

参考:
https://www.uniprot.org/
http://swxxx.alljournals.cn/html/2019/3/20190301.html#outline_anchor_15

UniprotKB：即Uniprot知识数据库。分为Swissprot（人工注释的高质量数据库）和TrEMBL（计算机自动注释）两个子库，包含了一系列的蛋白质注释信息。
UniParc：即Uniprot蛋白质归档数据库。是强大的非冗余蛋白质数据库，无论是否为同一物种，无论来自哪个数据库，只要序列相同，就会被合并为同一条目。
UniRef: 即Uniprot蛋白质参考数据库。包括UniRef100（使用来自UniprotKB和UniParc的序列聚类，长度不小于11个氨基酸的相同序列和片段生成，序列间Identity为100）、UniRef90（使用来自UniRef100的序列聚类生成，序列间至少有90的Identity，并且与每个簇最长的序列有80％以上的overlap）、UniRef50（使用来自UniRef50的序列聚类生成，序列间至少有50的Identity，并且与每个簇最长的序列有80％以上的overlap）。
Proteome：即蛋白质组数据库。由于物种的不同组织、发育阶段、环境条件的转录组均不同，并且翻译后的处理、修饰不同，导致了不同的最终产物蛋白。蛋白质组数据库就是针对目前已经完全测序的物种所收集的全部蛋白质集合。

高级检索 && 下载

输入(family:globin OR androglobin) organism:"homo sapiens"，来指定蛋白家族为为globin，物种名为homo sapiens的条目。你可以在后面加上AND reviewed:yes来指定只搜索Swiss-Prot的高质量结果，或者点击左边筛选栏选择。这次我们搜索到了一共165条globin家族的结果，并且有13条高质量序列。注意：Androglobin在Uniprot上属于Peptidase C2 family，因此需要额外检索，同样你的基因家族序列也要注意这点，这里的androglobin并不是算在family里面的，如果你想指定加上Peptidase C2 family，可以写成(family:globin OR family:"Peptidase C2 family") organism:"homo sapiens"，但这样会检索到额外的Peptidase C2 family成员。(或者你可以搜family:"Peptidase C2 family" AND androglobin)

fig12.png

我们点击Download，可以看到有canonical和canonical & isoform两种选择，那么这两种选择到底有什么不同呢？

fig13.png

参考：https://www.uniprot.org/help/canonical_and_isoforms
canonical sequence：为了减少冗余，UniProtKB/Swiss-Prot的策略是将同一个物种的同一基因编码的全部蛋白产物描述在一个条目内，即选取一个代表蛋白，筛选的标准至少满足以下条件之一：①所选的代表蛋白最流行。②代表蛋白与在其它物种中发现的直系同源序列相似度最高。③代表蛋白凭借其长度或氨基酸组成，它可以最清晰地描述结构域、异构体、多态性、翻译后修饰等。④在没有任何信息的情况下，选择最长的序列。
各种UniProtKB发布格式（Flat Text, XML, RDF/XML）只显示canonical序列。但是请注意！一些来自同一基因的替代性剪接亚型只共享几个外显子，如果有的话，对于一些'trans-splicing'事件也是一样的。在这些情况下，分歧显然太重要了，无法将所有的蛋白质序列合并到一个条目中，必须在单独的 "外部 "条目中描述异构体。由于来自大规模测序项目的序列数据增加，UniProtKB/TrEMBL可能包含额外的预测序列，这些序列由UniProtKB/Swiss-Prot条目中描述的基因编码。
根据上述的了解，以及经过本人多次验证，如果你是下Swissprot的序列，并且点上canonical的选项，那么只会下载canonical序列，而不会下isoform。但是如果你是要下trembl的序列，那么就没办法把单独一个基因对应的代表的唯一isoform下下来，也就是说，我们最好的做法还是尽可能选择模式物种的swissprot的序列作为参考，去下载。因为非模式物种的人工审核的序列实在太少，但是trembl序列质量又不高而且还有冗余。因此我们这里最好的做法就是选择swissprot并且点上canonical下载，即可得到需要的高质量人类globin序列集合