Monthly Archives: 6月 2021

六 24

什么才是生物信息学呢？

Posted on by

看到了朋友圈有人转发了这样的一个系列：

创作者的心声很有意思：想和大家一起开始学习生物信息学。因为最近看文献（涉及生信），总是大概明白了，自己还是无法操作。我也关注了很多生信的公众号，加了几个群，心动种草了一些线上课程。但是实在是觉得自己的基础知识太薄弱了，都看不懂群里的提问。所以决定要系统学一下，先看几本书吧，然后再学个线上操作课。

让我想起来了一个故事，就是现阶段生物信息学本科培养的人才与社会工作实践脱钩的问题。因为我看到他这个学习内容完全就是现阶段生物信息学本科培养大纲。

比如，他最新的一个是：【陪你学·生信】十二、RNA相关的简单分析，提到了Mfold 可以做RNA的二级结构预测，很有意思。在其网页工具上面输入你需要fold的序列，点击 fold RNA。如果事先不知道关于这段序列的任何信息，那么其他的参数都保持默认，就可以拿到该RNA序列的二级结构，确实是超级好用。

其实呢，这样的绝大部分生物信息学本科专业所教的知识点啦，但是这些呢，并不是社会实践中的生物信息学，真正的生物信息学不仅仅是能使用它这个mfold的网页工具，还需要能解决问题。比如这个网页工具是有序列输入长度限制的，如果你的序列超级长，比如最近流行的新冠病毒的某个RNA基因，就会超出限制，那么我们就得使用另外一个方法来进行mfold， 就是命令行软件使用。

软件安装3部曲

很久以前我分享过1000个生物信息学软件安装，见：生物信息学常见1000个软件的安装代码！，它们是有规律的，比如对C语言的源代码软件来说，就是3部曲啦，首先，登陆你的Linux系统，然后敲入下面的命令，全部的代码如下：

mkdir -p ~/biosoft/mfold && cd ~/biosoft/mfold 
wget http://www.unafold.org/download/mfold-3.6.tar.gz
tar zxvf mfold-3.6.tar.gz
cd mfold-3.6
./configure prefix=$HOME/biosoft/bin/
make
make install
# 前提是你的电脑存在 $HOME/biosoft/bin/ 文件夹，并且该 $HOME/biosoft/bin/ 文件夹是在环境变量
# 如果不在环境变量，就需要下面的2行代码哦！
PATH=$PATH:$HOME/biosoft/bin
export PATH

这样的话，你的Linux系统里面就有了mfold这个命令，如果你认真看它，其实就会发现它是一个shell脚本：

#!/bin/bash
# This shell script folds an RNA or DNA sequence and creates output
# files.
export _POSIX2_Version=0
export Package_URL="http://mfold.rna.albany.edu"

export DATDIR=`mfold_datdir|sed -e 's@/$@@'`
export Package=mfold
export Version=3.6

# Abort subroutine
abort() {
 rm -f mfold.log fort.*
 if [ $# -gt 0 ] ; then
 echo -e "$1"
 fi
 echo "Job Aborted"
 exit 1
}

虽然如此简单，但是没有这个基础的小伙伴仍然是会花式报错哦：

image-20210111214306533

输入数据的准备

既然我们在其网页工具都是输入序列即可，软件使用理论上也是如此，前面的网页工具使用的时候，我就是输入了软件安装包里面自带的测试序列：

>MDV-1 (-) RNA
GGGGAACCCCCCUUCGGGGGUCACCUCGCGCAGCGGGCUGCGCGAAGGGGCCACGCUGCGAAGCAGCGUG
GCGGUUCUCGUGCGUUACCGAAACGCACGAAGGUCGCGCCUCUUCACGAGGCGUCACCUGGGAGAGCGCG
AAAGCGCUAGCCCGUGACUCGUCACGGUCGAACUCCCGUACGAGGUGCCCGCACCUCGUCCCCCCUUCCG
GGGGGUCCCCA

默认参数运行mfold

如果前面的对mfold软件的C语言的源代码安装三部曲正常的话，该软件就可以调用，先看看帮助文档，如下：

$mfold
Usage is
mfold SEQ='file_name' with optional parameters:
 [ AUX='auxfile_name' ] [ RUN_TYPE=text (default) or html ]
 [ NA=RNA (default) or DNA ] [ LC=sequence type (default = linear) ]
 [ T=temperature (default = 37 deg C) ] [ P=percent (default = 5) ]
 [ NA_CONC=Na+ molar concentration (default = 1.0) ]
 [ MG_CONC=Mg++ molar concentration (default = 0.0) ]
 [ W=window parameter (default - set by sequence length) ]
 [ MAXBP=max base pair distance (default - no limit) ]
 [ MAX=maximum number of foldings to be computed (default 100) ]
 [ MAX_LP=maximum bulge/interior loop size (default 30) ]
 [ MAX_AS=maximum asymmetry of a bulge/interior loop (default 30) ]
 [ ANN=structure annotation type: none (default), p-num or ss-count ]
 [ MODE=structure display mode: auto (default), bases or lines ]
 [ LAB_FR=base numbering frequency ] [ ROT_ANG=structure rotation angle ]
 [ START=5' base # (default = 1)] [ STOP=3' base # (default = end) ]
 [ REUSE=NO/YES (default=NO) reuse existing .sav file ]

对于软件安装包里面自带的测试序列，我们可以这样运行：

$mfold SEQ='test.txt' 
# 不知道为什么我的报错了，如下：
mfold version 3.6
REUSE= NO
test.txt.pnt created.
Sequence length is 221
Folding at 37 degrees using version 3.0 dat files.
Save file created using nafold.
Save file is empty. No foldings.
Job Aborted

好奇怪啊，我也不知道为什么我的报错了。因为这个软件并不是我的刚需，我也懒得去解决这个bug啦，如果有这方面经验的小伙伴欢迎留言交流哈！

输出结果解读

如果是正确运行，大概会出现如下结果：

542K Mar 21 2009 mdv1.37.ct
5.1K Mar 21 2009 mdv1.37.ext
 13K Mar 21 2009 mdv1.37.plot
 39 Mar 21 2009 mdv1.dG
 61K Oct 24 2009 mdv1.jpg
 239 Nov 19 2009 mdv1-local.seq
 33 Nov 19 2009 mdv1.log
 51K Oct 24 2009 mdv1.pdf
 14K Mar 21 2009 mdv1.plot
 12K Oct 24 2009 mdv1.png
 21K Oct 24 2009 mdv1.ps
 109 Mar 21 2009 mdv1.run
4.0K Jan 11 16:12 mfold

相当于你在网页工具里面提交序列，选择参数，然后运行拿到的交互式结果。

其它编程语言的mfold

比如：https://pypi.org/project/seqfold/， seqfold is an implementation of the Zuker, 1981 dynamic programming algorithm, the basis for UNAFold/mfold, with energy functions from SantaLucia, 2004 (DNA) and Turner, 2009 (RNA).

类似的，也会有R语言，java等等开发的mfold。

conda或者docker

更高级啦，这里略。如果你感兴趣，可以考虑加入我们的《生信小成之conda交流群》。

我们邀请到了，简书conda教程单篇阅读量破40万的人气作者卖萌哥为咱们《生信技能树》和《生信菜鸟团》粉丝在钉钉群直播授课。直播是免费的哈，赶快下载钉钉软件加入吧，“Linux公益课(2021) 生物信息学”群的钉钉群号：33840083，下周六（2021-01-16）晚上八点开课哈。

同时我们提供一个微信交流群（钉钉软件我们并不是随时在线，不方便交流，钉钉仅仅是直播授课时候开启聊天），还是老规矩，18 元进群，一个简单的门槛，隔绝那些营销号！，仅此而已，考虑清楚哦！ 进群方式详见公众号推文：很多事情不一定有答案（但是可以有交流渠道）

我的疑问是

到底是各式各样的网页工具操作是生物信息学呢，还是基于Linux的NGS操作才是生物信息学呢？

​
六 24

为什么百分百还原文献结果反而不对呢

Posted on by

去年我在《生信技能树》给学员分享过一个案例：GSE16561 复现，小众的illumina bead芯片，链接在：http://www.bio-info-trainee.com/7483.html ，当时大费周折才拿到跟其数据集原始文献类似的结果，已经是很满意了。最近把这个数据集作为任务安排给最新学徒们，他们反馈给我的结果让我丈二和尚摸不着头脑，居然是百分百还原文献结果，如下所示的差异基因列表： Continue reading →

六 22

使用curatedTCGAData下载TCGA数据库信息好用吗

Posted on by

好久没有写TCGA数据库教程了，因为TCGA计划早在2017年就陆陆续续停止了，我那个时候写了几百个教程并且录制了视频。

Continue reading →

六 22

你还缺乳腺癌表达量数据集吗

Posted on by

生存分析你还是在TCGA吗？

最近有粉丝求助说他研究乳腺癌做了单细胞转录组数据，定位到了一个稀有细胞亚群，先看它感兴趣的亚群细胞特异性基因的临床意义，问我有没有除了TCGA数据库之外的其它数据库资源推荐。恰好我做这方面就顺手检索了一下，发现了 curatedBreastData 包，值得推荐！

Continue reading →

六 22

基因集的转录因子富集分析

Posted on by

一般来说，大家拿到了感兴趣的基因集后，通常是做超几何分布检验看看富集到了什么生物学功能数据库，比如KEGG或者GO数据库，或者走gsea/gsva这样的富集分析，也是注释生物学功能数据库。 大家读我的表达芯片的公共数据库挖掘系列推文应该是够多了：

Continue reading →