蛋白质分类相关资源(持续更新)

小妖儿

有关CD-HIT
用来对聚类生物序列减少序列冗余，可以用于大规模数据的序列聚类，利用其构建非冗余（序列）蛋白数据库以提高搜索效率。

（1）CD-HIT软件下载地址：http://www.bioinformatics.org/project/filelist.php?group_id=350
基于cd-hit的主算法还有几个附属的工具：
cd-hit-2d比较两个蛋白数据库并报告其中的相似比对；
cd-hit-est对DNA/RNA数据库进行聚类；
cd-hit-2d比较两个核酸数据库。在下载的时候就都包括在内了。但暂时没有用到。

（2）在linux下操作：
tar xvf cdhit.gz
cd cdhit
make

（3）编译完成后把需要操作的文件放在cdhit目录下（cd-hit命令前一定要加./)
./cd-hit -i a.fa -o output.fa -c 0.4  -n 2
其中a.fa是你要输入的文件，output.fa是你要输出的文件，-c为sequence identity threshold，根据要求或者目的选择，默认为0.9，选择了-c之后-n会随之改变，如下：

-n 5 for thresholds 0.7 ~ 1.0

-n 4 for thresholds 0.6 ~ 0.7

-n 3 for thresholds 0.5 ~ 0.6

-n 2 for thresholds 0.4 ~ 0.5  

还有一些其他参数，诸如-M 使用的内存大小，-T使用的进程数。

（4）会有3个文件生成：
.fa  
.fa.bak.clstr
.fa.clstr
只有当所有的序列都比对完之后 .fa文件中才会有内容，即为所需文件
相似度0.4以下的CD-HIT参见
http://bbs.malab.cn/forum.php?mo ... =1464&fromuid=3

花思谢

刘滨老师2014年的论文提供了两种提取氨基酸序列特征的方法。
他提供了一个web：web地址
如果不能正常使用的话，可以下载附件中的源码，里面有现成的程序可以直接调用。也有源码，可以自己修改。使用方法见文件里的readme。
我把其中一种改成java版的，在source code子文件夹下面，可做参考。

zouquan

极限学习

zouquan

菜鸟请看这个帖子：以亚细胞定位为例

蛋白质的亚细胞定位  （蛋白质是分子，我们搞信息的人看就是个氨基酸字符串，跟DNA类似。蛋白质在细胞里到处都有，细胞你记得吧？有细胞质、细胞核、高尔基体、叶绿体等等，现在给你一个蛋白质序列，请用软件预测它是位于细胞质里，还是细胞核里，还是高尔基体里？还是到处都有？） （这工作很重要，药就是蛋白质，如果细胞核里有问题了，你吃的药只能到达细胞质，那岂不是白吃了？）

问题明白了吧？

下面介绍方法，这类方法统称“机器学习”(machine learning)方法。（虽然叫机器学习，但是和硬件一点关系没有，完全是算法程序）

给你个蛋白质序列ACMNKGYW，鬼知道他在细胞核里，还是叶绿体中。怎么办呢？你得先给点已知的样例啊，看看已知的样例有没有啥规律，然后才好预测未知的，对不对？这些已知的样本就称为“训练集”(training set)，那些要预测的就称为“测试集”(testing set)。

好了，给你点“训练集”，比如：
高尔基体里有下面这2个蛋白质序列：(>开头那行是蛋白质的名字，下面的就是蛋白序列)
>Q6IBS0
MAHQTGIHATEELKEFFAKARAGSVRLIKVVIEDEQLVLGASQEPVGRWDQDYDRAVLPL
LDAQQPCYLLYRLDSQNAQGFEWLFLAWSPDNSPVRLKMLYAATRATVKKEFGGGHIKDE
LFGTVKDDLSFAGYQKHLSSCAAPAPLTSAERELQQIRINEVKTEISVESKHQTLQGLAF
PLQPEAQRALQQLKQKMVNYIQMKLDLERETIELVHTEPTDVAQLPSRVPRDAARYHFFL
YKHTHEGDPLESVVFIYSMPGYKCSIKERMLYSSCKSRLLDSVEQDFHLEIAKKIEIGDG
AELTAEFLYDEVHPKQHAFKQAFAKPKGPGGKRGHKRLIRGPGENGDDS
>Q9UJV3
METLESELTCPICLELFEDPLLLPCAHSLCFSCAHRILVSSCSSGESIEPITAFQCPTCR
YVISLNHRGLDGLKRNVTLQNIIDRFQKASVSGPNSPSESRRERTYRPTTAMSSERIACQ
FCEQDPPRDAVKTCITCEVSYCDRCLRATHPNKKPFTSHRLVEPVPDTHLRGITCLDHEN

细胞核膜里面有这3个蛋白质序列：
>O95866
MAVFLQLLPLLLSRAQGNPGASLDGRPGDRVNLSCGGVSHPIRWVWAPSFPACKGLSKGR
RPILWASSSGTPTVPPLQPFVGRLRSLDSGIRRLELLLSAGDSGTFFCKGRHEDESRTVL
HVLGDRTYCKAPGPTHGSVYPQLLIPLLGAGLVLGLGALGLVWWLHRRLPPQPIRPLPRF
APLVKTEPQRPVKEEEPKIPGDLDQEPSLLYADLDHLALSRPRRLSTADPADASTIYAVV
V
>Q7Z602
MPGHNTSRNSSCDPIVTPHLISLYFIVLIGGLVGVISILFLLVKMNTRSVTTMAVINLVV
VHSVFLLTVPFRLTYLIKKTWMFGLPFCKFVSAMLHIHMYLTFLFYVVILVTRYLIFFKC
KDKVEFYRKLHAVAASAGMWTLVIVIVVPLVVSRYGIHEEYNEEHCFKFHKELAYTYVKI
INYMIVIFVIAVAVILLVFQVFIIMLMVQKLRHSLLSHQEFWAQLKNLFFIGVILVCFLP
YQFFRIYYLNVVTHSNACNSKVAFYNEIFLSVTAISCYDLLLFVFGGSHWFKQKIIGLWN
CVLCR
>Q8N1M1
MTVTYSSKVANATFFGFHRLLLKWRGSIYKLLYREFIVFAVLYTAISLVYRLLLTGVQKR
YFEKLSIYCDRYAEQIPVTFVLGFYVTLVVNRWWNQFVNLPWPDRLMFLISSSVHGSDEH
GRLLRRTLMRYVNLTSLLIFRSVSTAVYKRFPTMDHVVEAGFMTTDERKLFNHLKSPHLK
YWVPFIWFGNLATKARNEGRIRDSVDLQSLMTEMNRYRSWCSLLFGYDWVGIPLVYTQVV

其他的位置也给了类似的样例。好了，来个测试样例：
>P60880
MAEDADMRNELEEMQRRADQLADESLESTRRMLQLVEESKDAGIRTLVMLDEQGEQLERI
EEGMDQINKDMKEAEKNLTDLGKFCGLCVCPCNKLKSSDAYKKAWGNNQDGVVASQPARV
VDEREQMAISGGFIRRVTNDARENEMDENLEQVSGIIGNLRHMALDMGNEIDTQNRQIDR
IMEKADSNKTRIDEANQRATKMLGSG

请问他位于哪个位置？

你用肉眼是看不出来，对不对？跟谁也不太像，那怎么办？字符串比较像不像太难了。

“机器学习”的核心思想就是把这些难于比较的东西（比如字符串、语音信号、图片等等），统统转化为向量（向量就是矩阵的一行，就是几个数）。

字符串比较难计算相似度，但向量很容易啊，欧氏距离就可以了，对不对？

把字符串（或者图片、语音信号）转化成向量的过程就称为“特征提取”(feature extraction)，这个向量就称为特征(feature)

怎么提啊？最简单的：统计一下每个字母出现的频率。蛋白质序列中只有20个字符（注意不是26个，氨基酸只有20种）每个字母出现的频率计算一下（比如A出现了10次，序列长度为100，那就是0.1咯），一条序列就变成了一个20维的向量（就是20个数值），对不对？

你把训练集中所有的样本都转化成20维向量，测试集的样本也转化成向量，计算一下他们的欧氏距离，测试集的样本和训练集中的哪个最近，就猜是跟他一样的位置。可以吧？这种预测（分类）方法就称为 “最近邻”。

当然这种分类算法比较土，还有很多高大上的算法，如：支持向量机、随机森林、集成学习。


好了，如果明白了，自己上网查一查以下概念：
分类、训练集、测试集、样本、特征（属性）、交叉验证（这个我没讲，自己理解一下）


后续工作：

1. 数据收集。要整理好细胞一共有多少个位置，每个位置上都有哪些蛋白质，有没有哪些蛋白质出现在多个位置？

2. 算频率这个特征提取算法太土了，有很多高大上的办法，前面的帖子有很多现成的源代码，只需要运行就可以。

3. 最近邻这种分类器也太土了，也有不少算法，也有现成的软件，如weka。

4. 你的工作就是尝试，哪种特征配合哪种分类器效果最好？试出来就行。这就叫做实验。




多说几句：

别觉得蛋白质这玩意没意思又没用，实际上在教你如何用“机器学习”的思想进行预测，“黄色图片识别”、“语音识别”、“指纹识别”、“复杂仪器的故障诊断”用的都是一样的方法。你学会了这个，可以去干你感兴趣的事情。只不过那些工作要么数据难获得、要么论文不好发。菜鸟阶段，先干点经济实惠的，学有所成了你再去整高大上的。

RockRabbit

MEME软件可以帮助你找数据集中的functional motifs或者conservation motifs，可以通过这些motifs寻找些特征。供大家学习。

http://meme-suite.org/tools/meme

zouquan

二类分类（特殊蛋白识别）与PSI-BLAST比较效果

取出CD-HIT之后的正反例的fasta文件。可以分别修改一下每条序列的名字，使得从名字就能一下看出是正反例，然后把正反例合并在一起。

    用合并在一起的fasta文件当成数据库，每一条序列依次当成查询文件，用PSI-BLAST进行查询，肯定自己和自己得分最高，得分第二高的如果和查询序列类别相同（同是正例，或同是反例），则正确，反之则错误。统计一下正反例的混淆矩阵。

RockRabbit

一个在线的蛋白质各种特征提取网站：

http://bcb.ncat.edu/Features/

使用注意：输入的文件名必须得是.fasta结尾的。如果是做二类分类问题，正例文件可以命名为group1.fasta，反例文件可以命名为group2.fasta，所提取到的特征里的最后一项为他默认给你添加的类标签，例如group1.fasta所对应的特征文件中的类标签为0，group2.fasta所对应的特征文件中的类标签为1。多分类问题依次类推。

RockRabbit

能生成21种基于PSSM的特征：

http://possum.erc.monash.edu/

zouquan

188D代码具体特征意义：
(有了这个解释，可以用降维的方法，降到最低，从而看看到底哪些特征可以对该蛋白做最有效的分类)

前20维，分别是20种氨基酸（按字母序ACDEFGHIKLMNPQRSTVWY），在序列中的含量。（出现个数/序列长度）

21-41维是疏水性特征
        /**   calc_Hydrophobic();
         * string sP="RKEDQN";亲水
         * string sN="GASTPHY";中性 （H应该去亲水）
         * string sH="CVLIMFW";疏水  （W应该去中性）
         */

21-23维分别是亲水、中性、疏水的氨基酸含量 （出现个数/序列长度）

24-26维分别是转换频率，亲水/中性，亲水/疏水，中性/疏水，（出现转换次数/序列长度-1）

27-31维分别是亲水氨基酸第1个，25%个，50%个，75%个和最后一个在序列中的位置，（第n位/序列长度）
32-36维分别是中性氨基酸第1个，25%个，50%个，75%个和最后一个在序列中的位置，（第n位/序列长度）
37-41维分别是疏水氨基酸第1个，25%个，50%个，75%个和最后一个在序列中的位置，（第n位/序列长度）

下面特征类似

42-62是范德华力
/**  calc_Vanderwaal()
         * string sP="GASCTPD";
         * string sN="NVEQIL";
         * string sH="MHKFRYW";
         */

63-83极性
        /**  calc_Polarity1()
         * string sP="LIFWCMVY";
         * string sN="PATGS";
         * string sH="HQRKNED";
         */

84-104是极化性质
        /**  calc_Polarizability1()
         * string sP="GASDT";
         * string sN="CPNVEQIL";
         * string sH="KMHFRYW";
         */

105-125是电荷性质
        /** calc_Charge()
         * string sP="KR";
         * string sN="ANCQGHILMFPSTWYV";
         * string sH="DE";
         */

126-146是表面张力
        /**  calc_Surfacetension()
         * string sP="GQDNAHR";
         * string sN="KTSEC";
         * string sH="ILMFPWYV";
         */

147-167是二级结构
        /**  calc_Secondarystructure()
         * string sP="EALMQKRH";
         * string sN="VIYCWFT";
         * string sH="GNPSD";
         */
168-188是溶剂可及性
        /**  calc_Solventaccessibility()
         * string sP="ALFCGIVW";
         * string sN="RKQEND";
         * string sH="MPSTHY";
         */