0%

写在前面

长非编码RNA(Long non-coding RNA, lncRNA)是长度大于 200 个核苷酸的非编码 RNA。200个核苷酸的限制可以区分当时已知的功能非编码RNA(如miRNA、miRNA前体)和结构非编码RNA(如tRNA、rRNA)。可想而知,lncRNA是个大杂烩,种类繁多,功能复杂。

随着转录组测序的大规模应用,大量的lncRNA被发现。如果你做了一个转录组研究,可能也会有鉴定新lncRNA的需求。转录组构建之后,去掉已知的编码蛋白基因和lncRNA,留下的就是新发现的转录本。从新发现的转录本鉴定lncRNA,简单来讲,就是去掉转录噪音和去掉编码蛋白基因。

转录组测序可以分析出几万到几十万转录本,其中很大一部分是转录噪音。我们通常来讲,会去掉单外显子和低表达转录本。有人认为,RNA聚合酶II可能随机结合在DNA上,瞬时转录出一些RNA,因为不是正常的基因转录,所以会马上关闭。这些转录噪音的特征就是单外显子和低表达。

通常根据lncRNA在基因组上与蛋白质编码基因的相对位置将其分为4类:基因间长链非编码RNA(intergenic lncRNA,由基因组上编码基因之间的间隔序列转录产生,也即long intergenic ncRNA, lincRNA)、内含子长链非编码RNA(intronic lncRNA,完全由蛋白质编码基因内的内含子转录而来)、正义链非编码RNA(sense lncRNA,转录方向与蛋白质编码基因转录方向相同,并且与蛋白质编码基因外显子有部分或者全部重叠)、反义长链非编码RNA(antisense lncRNA,转录方向与蛋白质编码基因转录方向相反)。其中lincRNA争议很小,sense lncRNA和intronic lncRNA可能是基因转录错误或者转录加工得到中间产物。如果是链非特异性的库,转录本的方向是软件推断出来的,antisense lncRNA可能就是sense lncRNA,也可能是转录噪音。有些严格的一点的研究,可能只要intergenic的新转录本。

接下来是去掉编码蛋白基因,因为lncRNA特征不明显,鉴定lncRNA就变成鉴定编码蛋白基因,然后取反。有人自己编写程序去掉编码蛋白基因,比如BLAST比对已知蛋白库(不能比对上或者限制比对上的长度),再加上开放阅读框长度不能大于30(或者某个其他值)。鉴定lncRNA的工具其实很多,大家常用的套路也是综合多个工具的结果。本文推荐CPC2和CPC这两个软件。

使用CPC2鉴定lncRNA

CPC2只使用序列本身的特征,速度非常快。我现在的做法是先跑CPC2,去掉coding gene,剩下的再去跑CPC。CPC2的官方地址:http://cpc2.cbi.pku.edu.cn/download.php

a. CPC2的安装

CPC2依赖Biopython,所以我们先安装它。

1
pip install biopython

然后安装CPC2。

1
2
3
4
5
6
7
tar zxvf CPC2-beta.tar.gz
cd CPC2-beta
export CPC_HOME="$PWD"
cd libs/libsvm
tar zxvf libsvm-3.18.tar.gz
cd libsvm-3.18/
make clean && make

注意CPC2目前只支持Python2,如果你装了Anaconda,可以开一个Python2的虚拟环境。export CPC_HOME="$PWD"这句是将CPC_HOME这个环境变量设置为CPC2的主目录路径,重启终端后,这个设置就无效了。如果要永久设置可以修改~/.bashrc文件,添加export CPC_HOME="/home/chenwen/bin/CPC2-beta",注意你的实际路径可能和我的不同。

b. 运行CPC2

1
2
cd $CPC_HOME
python ./bin/CPC2.py -i input.fa -o output.txt

使用CPC鉴定lncRNA

CPC是鉴定lncRNA经典软件,它先把未知转录本和已知蛋白库BLAST对比,BLAST的结果作为libsvm(一个机器学习分类库,CPC2也用了它)的分类特征。因为有了BLAST,可靠性更高,同时也非常慢。

CPC是个很老的软件了,它调用的是blastall,也就是老版本的BLAST,而不是新版本的BLAST+。CPC和老版本BLAST的下载地址,有些朋友可能找不到,我把它们直接放在我的网站上了。

CPC下载地址:https://static.biochen.com/software/cpc-0.9-r2.tar.gz
老版本的BLAST下载地址:https://static.biochen.com/software/blast-2.2.26-x64-linux.tar.gz

a. CPC的安装

1
2
3
4
5
6
7
8
9
tar -zxvf cpc-0.9-r2.tar.gz
cd cpc-0.9-r2/libs/libsvm
tar -zxvf libsvm-2.81.tar.gz
cd libsvm-2.81
make clean && make
cd ../..
tar -zxvf estate.tar.gz
cd estate
make clean && make

b. 建立BLAST库

需要使用蛋白质库,UniRef90或者NCBI的nr都可以,用formatdb命令建库时,必须命名为”prot_db”, 且放在CPC安装目录下的data目录下面。

1
2
cd cpc-0.9-r2/data
formatdb -i (your_fasta_file) -p T -n prot_db

c. 运行CPC

1
2
cd cpc-0.9-r2/
bin/run_predict.sh (input_seq) (result_in_table) (working_dir) (result_evidence)

括弧里面的内容需要改成你实际的文件名或者路径。run_predict.sh会调用远程blast,建议运行run_predict_local.sh,并把这个文件中blast_opts=”$blast_opts -a 2″; # 2CPUs, boost the performance这句话中的2,改成你实际电脑使用的CPU核数。

有些命令用得多了,也就记住了,比如cd和ls。有些命令用得没那么多,记不住,拿个小本本记下来,于是有了这篇博客。

压缩和解压

.tar.gz压缩

1
tar -zcvf compressed.tar.gz to_compress_file_or_folder

-p 保留文件权限

.tar.gz解压

1
tar -zxvf compressed.tar.gz

.zip压缩

1
2
zip -r compressed.zip to_compress_folder
zip compressed.zip to_compress_file

.zip 解压

1
unzip compressed.zip

.gz压缩

1
gzip to_compress_file

.gz解压

1
gunzip compressed.gz

.tar.bz2解压

1
tar -jxvf compressed.tar.bz2

.tar.xz解压

1
tar -xvJf compressed.tar.xz

之所以没有压缩为.tar.bz2和.tar.xz的命令,是因为这不是命令大全,只记录我用过的命令吧。Linux文件压缩首选.tar.gz格式,加-p参数可以保留文件权限信息。如果为了考虑Windows也能很方便打开,可以选择压缩为.zip格式。

添加环境变量

1
2
echo 'export PATH=/home/biochen/bin:$PATH' >> ~/.bashrc
source ~/.bashrc

很多Linux软件提供了编译好的二进制文件,可以在软件的目录输入./software_name来使用它。如果要想到处可以用,将软件的路劲添加到环境变量即可。

scp

上传

1
2
scp -P port local_file remote_username@remote_ip:remote_folder
scp -P port -r local_folder remote_username@remote_ip:remote_folder

下载

1
2
scp -P port remote_username@remote_ip:remote_file local_folder
scp -P port -r remote_username@remote_ip:remote_folder local_folder

自己的电脑与远程服务器的上传下载,我一般使用FileZilla,鼠标点几下还是很香的。如果要在两台服务器之间传数据,特别是大一点的数据,下载到自己的电脑,然后上传,就没有那么香了,这时候scp命令就可以派上用场了。

查看进程

1
2
ps -elf |grep user_name
ps -elf |grep software_name

查看进程命令是ps -elf,会输出一堆进程,后面一般通过管道接grep,筛选出特定用户或者特定软件的进程。进程没什么好看的哈,主要是为了找出进程的id(第二列的数字),然后用kill命令杀掉。

1
kill job_id

查看网络端口

1
sudo netstat -tunlp |grep 80

或者

1
sudo lsof -i:80

以上命令,找出哪个进程监听80端口,然后等待这个进程的命运一般也是kill。

修改文件所有者和权限

修改文件的所有者

1
2
chown <user> <file|directory> [-R]
chgrp <group> <file|directory> [-R]

尖括弧<>括起来的参数是必须要的,方括弧[]括起来的参数是可选的,竖杠|两边二选一。-R递归文件夹下所有文件和文件夹,不然只是修改文件夹本身的所有者。

修改文件的权限

1
chmod 777 <file|directory> [-R]

Linux的有着比Windows复杂的文件权限,然而因为权限设置不对会产生很多问题,无脑设置为777(全部权限)是我通常的做法。Linux的文件权限机制的初衷是为了安全,如果我知道该给多少权限时,还会改过来。以下是常用的权限码:

1
2
3
4
5
6
7
-rw------- (600)  只有拥有者有读写权限
-rw-r--r-- (644) 只有拥有者有读写权限;而属组用户和其他用户只有读权限
-rwx------ (700) 只有拥有者有读、写、执行权限
-rwxr-xr-x (755) 拥有者有读、写、执行权限;而属组用户和其他用户只有读、执行权限
-rwx--x--x (711) 拥有者有读、写、执行权限;而属组用户和其他用户只有执行权限。
-rw-rw-rw- (666) 所有用户都有文件读、写权限。
-rwxrwxrwx (777) 所有用户都有读、写、执行权限。

查看文件

查看文件前10行

1
top -n 10 file_name

查看文件后10行

1
tail -n 10 file_name

查看文件11-20行

1
sed -n '11,20p file_name

查看文件第21到最后一行

1
sed -n '21,$p' file_name

删掉文件前10行

1
sed -i '1, 10d' file_name

查看文件有多少行

1
wc -l file_name

查看文件有多少列

1
awk -F '\t' 'NR==1 {print NF}' file_name

查看表头

1
less file_name |head -n 1 |tr '\t' '\n' |cat -n

未完待续…

从NCBI下载对应系统的BLAST+程序

根据自己的操作系统,从NCBI的FTP下载合适的BLAST+版本。我的操作系统是Ubuntu 18.04,选择的是ncbi-blast-2.10.0+-x64-linux.tar.gz。

1
ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/

解压,将程序路径添加到环境变量

1
2
tar –zxvf ncbi-blast-2.10.0+-x64-linux.tar.gz
vi ~/.bashrc

在末尾添加:

1
export PATH=/home/chenwen/bin/ncbi-blast-2.10.0+-x64-linux/bin:$PATH

这里的路径视具体情况而定,我喜欢将软件都放在/home/chenwen/bin目录下
更新,使配置生效:

1
source ~/.bashrc

从NCBI下载数据库

下载地址:

1
ftp://ftp.ncbi.nih.gov/blast/db/

有关每个文件的含义请阅读README文件。下面摘录几种:
human_genomic.gz 人类基因组序列
nr.gz All non-redundant GenBank CDS translations + RefSeq Proteins + PDB + SwissProt + PIR + PRF所有非冗余的的GenBank CDS区的翻译序列 + 参考序列的蛋白 + PDB数据库 + SwissProt蛋白数据库 + PRF蛋白数据库
nt.gz 除wgs, gss, sts, pat, est, htg以外的核酸序列,注意不是非冗余的
htg.gz 来源于GenBank, EMBL, and DDBJ的高通量基因组测序序列

格式化数据库

BLAST+使用makeblastdb命令格式化数据库。老版本BLAST使用formatdb命令。

1
makeblastdb -in db.fasta -dbtype prot -parse_seqids -out db_name

参数说明:
-in:待格式化的序列文件
-dbtype:数据库类型,prot或nucl
-out:数据库名

运行BLAST+

BLAST+提供多种比对程序:
blastp: 用蛋白质序列搜索蛋白质序列库
balstn: 用核酸序列搜索核酸库
blastx:核酸序列对蛋白库的比对,核酸序列在比对之前自动按照六个读码框翻译成蛋白质序列
tblastn:蛋白质序列对核酸库的比对,核酸库中的序列按照六个读码框翻译后与蛋白质序列进行比对搜索
tblastx:核酸序列对核酸库在蛋白质级别的比对,两者都在搜索之前翻译城蛋白质进行比对
核酸序列比对核酸数据库(blastn):

1
blastn -query seq.fasta -out out.txt -outfmt 6 -db dbname -evalue 1e-5 -num_threads 8

参数说明:
-query: 输入文件路径及文件名
-out:输出文件路径及文件名
-db:格式化了的数据库路径及数据库名
-evalue:设置输出结果的e-value值
-num_threads:线程数
-outfmt:输出文件格式,6是tabular格式
-max_target_seqs:设置最大目标序列匹配数
其他程序比对跟blastn相似,更多参数可以用-help查询。

-outfmt 6的输出文件格式详解

1
2
3
4
5
6
7
8
9
10
11
12
[00] Query id
[01] Subject id
[02] % identity
[03] alignment length
[04] mismatches
[05] gap openings
[06] q. start
[07] q. end
[08] s. start
[09] s. end
[10] e-value
[11] bit score

安装MySQL

在Ubuntu中使用apt安装MySQL Server,如果是CentOS,将apt替换为yum。

1
sudo apt install mysql-server

以前安装过程中需要设置MySQL root用户密码,现在Ubuntu 18.04和CentOS 8安装MySQL都不需要设置root的密码了。安装完MySQL后,运行mysql_secure_installation,设置root密码和修改一些安全配置。

1
sudo mysql_secure_installation

启动、重启、查看、登录MySQL

启动MySQL

1
systemctl start mysqld

重启MySQL

1
systemctl restart mysqld

查看MySQL状态

1
systemctl status mysql.service

登录MySQL

1
mysql -u root -p

MySQL基础

MySQL语句约定使用大写,其实它是大小写不敏感的,语句以分号;结尾。
查看MySQL版本

1
SELECT VERSION();
1
2
3
4
5
6
+-----------+
| VERSION() |
+-----------+
| 8.0.17 |
+-----------+
1 row in set (0.01 sec)

查看数据库

1
SHOW DATABASES;
1
2
3
4
5
6
7
8
9
+--------------------+
| Database |
+--------------------+
| information_schema |
| mysql |
| performance_schema |
| sys |
+--------------------+
4 rows in set (0.01 sec)

SQL初步

新建数据库

1
CREATE DATABASE eg_db;

删除数据库

1
DROP DATABASE eg_db;

选定数据库
对表进行操作之前,需要选定数据库

1
USE eg_db;

创建表格
数据以表格形式保存,每一行是一个记录,每一列是一个属性。

1
2
3
4
5
CREATE TABLE person (
id int,
name varchar(50),
age int
);

查看表结构

1
DESCRIBE person;
1
2
3
4
5
6
7
8
+-------+-------------+------+-----+---------+-------+
| Field | Type | Null | Key | Default | Extra |
+-------+-------------+------+-----+---------+-------+
| id | int(11) | YES | | NULL | |
| name | varchar(50) | YES | | NULL | |
| age | int(11) | YES | | NULL | |
+-------+-------------+------+-----+---------+-------+
3 rows in set (0.00 sec)

删除表

1
DROP TABLE person;

重命名表

1
ALTER TABLE person RENAME people;

向表中添加一列

1
ALTER TABLE person ADD address varchar(100);

删除表中一列

1
ALTER TABLE person DROP COLUMN address;

修改一个列的数据类型

1
ALTER TABLE person MODIFY address varchar(50);

重命名一个列

1
ALTER TABLE person CHANGE COLUMN address addr varchar(50);

向表格中插入一条记录
INSERT INTO 表名称 VALUES (值1,值2,…)

1
INSERT INTO person VALUES (1, 'chenwen', 28, 'HNU');

或:INSERT INTO 表名称 (列1,列2,…) VALUES (值1,值2,…)

从表格中查询记录
SELECT 列名称1,列名称2… FROM 表名称;

1
SELECT name FROM person;

或:SELECT * FROM 表名称;

从表格中按条件查询一条记录
SELECT 列名称 FROM 表名 WHERE 列 运算符 值;

1
SELECT * FROM person WHERE name='chenwen';

从表格中删除一条记录
DELETE FROM 表名称 WHERE 列 运算符 值;

1
DELETE FROM person WHERE name='chenwen';

或:DELETE * FROM 表名称;

从表格中更新一条记录
UPDATE 表名称 SET 列名称=新值 WHERE 列=值;

1
UPDATE person SET age=31 WHERE name='chenwen';

返回结果删除重复项
SELECT DISTINCT 列名称 FROM 表名称;

1
SELECT DISTINCT name FROM person;

对查询结果按指定列进行排序
SELECT * FROM 表名称 ORDER BY 列名称;

1
SELECT * FROM person ORDER BY id;

或:SELECT * FROM 表名称 ORDER BY 列名称 DESC;

1
SELECT * FROM person ORDER BY id DESC;

以上命令很多,掌握它们的要领是按照增、删、改、查去运用它们,这也是我们实际上最常用的操作。操作数据库之前,记得使用USE选定这个数据库。

MySQL用户及权限管理

MySQL数据库默认只有一个root用户,MySQL将用户信息保存在mysql数据库user表中。
创建一个新用户
CREATE USER 用户名 IDENTIFIED BY ‘密码’;

1
CREATE USER chenwen IDENTIFIED BY 'password';

新用户不能创建新的数据库,因为没有设置权限。我们使用root账号创建新用户后,也顺便创建这个用户操作的数据库,并授予新用户权限,而不是用新用户去创建它的数据库。

删除用户
DROP USER 用户名;

1
DROP USER chenwen;

重命名用户
RENAME USER 原用户名 TO 新用户名;

1
RENAME USER chenwen TO cw;

修改当前用户密码
SET PASSWORD=PASSWORD(‘新密码’);

1
SET PASSWORD=PASSWORD('new_password');

修改指定用户密码
SET PASSWORD FOR 用户名 = PASSWORD(‘新密码’);

1
SET PASSWORD FOR chenwen = PASSWORD('new_password');

授予权限
授予一个用户权限:GRANT ALL PRIVILEGES ON 层级 TO 用户名@主机 WITH GRANT OPTION;
例如,授予chenwen用户针对eg_db数据库全部权限

1
GRANT ALL PRIVILEGES ON eg_db.* TO chenwen@'%' WITH GRANT OPTION;

如果你用的是老版本的MySQL,比如5.7,授权的命令不一样

1
GRANT ALL PRIVILEGES ON eg_db.* TO chenwen@'%' IDENTIFIED BY 'password';

可以使用以下方式指定权限范围

1
2
3
4
1. 所有主机:'%'
2. 精确的主机名或IP地址:www.biochen.com或者192.168.1.1
3. 使用"*"通配符:*.biochen.com
4. 指定一个网段:192.168.1.0/255.255.255.0

MySQL备份与恢复

备份数据库
mysqldump -u root -p 数据库名称 > 备份文件.sql

1
mysqldump -u root -p eg_db > eg_db.sql

mysqldump备份出来的是纯文本的SQL文件,可以修改后作为其他数据库数据使用。

从备份文件中恢复数据库
mysql -u root -p 数据库名称 < 备份文件.sql

1
mysqldump -u root -p eg_db > eg_db.sql

MySQL的字符编码

MySQL 8.0数据库的默认编码是utf8,老版本默认是latin,再也不用担心字符编码设置不对的乱码了。
以下命令可以查看MySQL当前使用的编码

1
SHOW VARIABLES LIKE 'character_set%';
1
2
3
4
5
6
7
8
9
10
11
12
13
+--------------------------+----------------------------+
| Variable_name | Value |
+--------------------------+----------------------------+
| character_set_client | utf8mb4 |
| character_set_connection | utf8mb4 |
| character_set_database | utf8mb4 |
| character_set_filesystem | binary |
| character_set_results | utf8mb4 |
| character_set_server | utf8mb4 |
| character_set_system | utf8 |
| character_sets_dir | /usr/share/mysql/charsets/ |
+--------------------------+----------------------------+
8 rows in set (0.00 sec)
1
SHOW VARIABLES LIKE 'collation%';
1
2
3
4
5
6
7
8
+----------------------+--------------------+
| Variable_name | Value |
+----------------------+--------------------+
| collation_connection | utf8mb4_0900_ai_ci |
| collation_database | utf8mb4_0900_ai_ci |
| collation_server | utf8mb4_0900_ai_ci |
+----------------------+--------------------+
3 rows in set (0.00 sec)

创建数据库时指定编码

1
2
3
CREATE DATABASE eg_db
DEFAULT CHARACTER SET utf8
DEFAULT COLLATE stf8_general_ci;

修改一个数据库的编码

1
ALTER DATABASE eg_db CHARACTER SET utf8 COLLATE utf8_general_ci;

MySQL 数据类型

MySQL中定义数据字段的类型对你数据库的优化是非常重要的。MySQL支持多种类型,大致可以分为三类:数值、日期/时间和字符串(字符)类型。
数值类型
MySQL支持所有标准SQL数值数据类型。这些类型包括严格数值数据类型(INTEGER、SMALLINT、DECIMAL和NUMERIC),以及近似数值数据类型(FLOAT、REAL和DOUBLE PRECISION)。关键字INT是INTEGER的同义词,关键字DEC是DECIMAL的同义词。作为SQL标准的扩展,MySQL也支持整数类型TINYINT、MEDIUMINT和BIGINT。下面的表显示了需要的每个整数类型的存储和范围。

类型 大小 范围(有符号) 范围(无符号) 用途
TINYINT 1 字节 (-128,127) (0,255) 小整数值
SMALLINT 2 字节 (-32 768,32 767) (0,65 535) 大整数值
MEDIUMINT 3 字节 (-8 388 608,8 388 607) (0,16 777 215) 大整数值
INT或INTEGER 4 字节 (-2 147 483 648,2 147 483 647) (0,4 294 967 295) 大整数值
BIGINT 8 字节 (-9 233 372 036 854 775 808,9 223 372 036 854 775 807) (0,18 446 744 073 709 551 615) 极大整数值
FLOAT 4 字节 (-3.402 823 466 E+38,1.175 494 351 E-38),0,(1.175 494 351 E-38,3.402 823 466 351 E+38) 0,(1.175 494 351 E-38,3.402 823 466 E+38) 单精度 浮点数值
DOUBLE 8 字节 (1.797 693 134 862 315 7 E+308,2.225 073 858 507 201 4 E-308),0,(2.225 073 858 507 201 4 E-308,1.797 693 134 862 315 7 E+308) 0,(2.225 073 858 507 201 4 E-308,1.797 693 134 862 315 7 E+308) 双精度 浮点数值
DECIMAL DECIMAL(M,D) ,如果M>D,为M+2,否则为D+2 依赖于M和D的值 依赖于M和D的值 小数值

日期和时间类型
表示时间值的日期和时间类型为DATETIME、DATE、TIMESTAMP、TIME和YEAR。每个时间类型有一个有效值范围和一个”零”值,当指定不合法的MySQL不能表示的值时使用”零”值。TIMESTAMP类型有专有的自动更新特性。

类型 大小 (字节) 范围 格式 用途
DATE 3 1000-01-01/9999-12-31 YYYY-MM-DD 日期值
TIME 3 ‘-838:59:59’/‘838:59:59’ HH:MM:SS 时间值或持续时间
YEAR 1 1901/2155 YYYY 年份值
DATETIME 8 1000-01-01 00:00:00/9999-12-31 23:59:59 YYYY-MM-DD HH:MM:SS 混合日期和时间值
TIMESTAMP 8 1970-01-01 00:00:00/2037 年某时 YYYYMMDD HHMMSS 混合日期和时间值,时间戳

字符串类型
字符串类型指CHAR、VARCHAR、BINARY、VARBINARY、BLOB、TEXT、ENUM和SET。CHAR和VARCHAR类型类似,但它们保存和检索的方式不同。它们的最大长度和是否尾部空格被保留等方面也不同。在存储或检索过程中不进行大小写转换。BINARY和VARBINARY类类似于CHAR和VARCHAR,不同的是它们包含二进制字符串而不要非二进制字符串。也就是说,它们包含字节字符串而不是字符字符串。这说明它们没有字符集,并且排序和比较基于列值字节的数值值。BLOB是一个二进制大对象,可以容纳可变数量的数据。有4种BLOB类型:TINYBLOB、BLOB、MEDIUMBLOB和LONGBLOB。它们只是可容纳值的最大长度不同。有4种TEXT类型:TINYTEXT、TEXT、MEDIUMTEXT和LONGTEXT。这些对应4种BLOB类型,有相同的最大长度和存储需求。

类型 大小 用途
CHAR 0-255字节 定长字符串
VARCHAR 0-65535 字节 变长字符串
TINYBLOB 0-255字节 不超过 255 个字符的二进制字符串
TINYTEXT 0-255字节 短文本字符串
BLOB 0-65 535字节 二进制形式的长文本数据
TEXT 0-65 535字节 长文本数据
MEDIUMBLOB 0-16 777 215字节 二进制形式的中等长度文本数据
MEDIUMTEXT 0-16 777 215字节 中等长度文本数据
LONGBLOB 0-4 294 967 295字节 二进制形式的极大文本数据
LONGTEXT 0-4 294 967 295字节 极大文本数据

结束

我大部分时候,只是创建一个MySQL数据库给其他程序调用,比如WordPress。以上的内容,基本适合这种轻量级应用了。感谢你看到结尾,最后我们以一个命令结束。
退出MySQL

1
EXIT;

对于实验生物学研究者来说,经常遇到的一个场景是,对照组测量3个值,实验组测量3个值,然后用T检验得出均值是否有显著性差异,然后做出结论。如下图所示:

那么T检验靠谱吗?本文采用模拟抽样的方法来探讨这件事情,由于本人统计学知识有限,还请各位读者批评指正。

我们通过构建两个随机总体,总体1的均值固定为1(模拟control),总体2的均值我们从依次取0.5到2.0(模拟treatment相对与control的变化),标准差我们取0.01、0.05、0.10(模拟测量误差1%,5%,10%)。总体数量为1000,每次随机抽取3个样本,抽样1000次,考察通过T检验得出正确结论的频率。代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
import numpy as np
import pandas as pd
import seaborn as sns
from scipy import stats
import matplotlib.pyplot as plt
plt.xkcd()
sns.set(style="white", context="talk")
def norm_t(loc1, loc2, scale1, scale2, n):
normSpace1 = stats.norm.rvs(loc = loc1, scale = scale1, size = 1000)
normSpace2 = stats.norm.rvs(loc = loc2, scale = scale2, size = 1000)
count = 0
for i in range(1000):
sample1 = np.random.choice(normSpace1, n)
sample2 = np.random.choice(normSpace2, n)
t, p = stats.ttest_ind(sample1, sample2)
if (p &lt; 0.05 and (sample1.mean() &gt; sample2.mean()) == (loc1 &gt; loc2)):
count = count + 1
return count / 1000
def plot(stds):
df = pd.DataFrame()
for std in stds:
x = np.arange(0.5, 2.1, 0.1)
y = [norm_t(loc1 = 1, loc2 = i, scale1 = std, scale2 = std, n = 3) for i in x]
z = [std for i in x]
df_std = pd.DataFrame([x, y, z]).T
df_std.columns = ["difference", "sensitivity", "std"]
if df.shape == (0, 0):
df = df_std
else:
df = df.append(df_std)
sns.pointplot(x = "difference", y = "sensitivity", hue = "std", data = df)
plot([0.01, 0.05, 0.10])

通过上图,我们发现,当标准差我们取0.01(模拟测量误差1%)时,增大到1.1倍或者减小到0.9倍,正确率可以达到100%。而当标准差取0.05和0.10(模拟测量误差5%和10%),效果就不理想了,本来是有差异的,还是有很大几率通过T检验得出错误的结论。显而易见的是,随着标准差的增加,T检验的效果变差了。

那么当标准差取0.10时,随着抽样样本数的增加,T检验判断正确的频率会是怎么样的呢?请看如下代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
def plot2(diffs):
df = pd.DataFrame()
for diff in diffs:
x = list(range(3, 60, 3))
y = [norm_t(loc1 = 1, loc2 =1.1, scale1 = 0.10, scale2 = 0.10, n = i) for i in x]
z = [diff for i in x]
df_diff = pd.DataFrame([x, y, z]).T
df_diff.columns = ["samle number", "sensitivity", "difference"]
if df.shape == (0, 0):
df = df_diff
else:
df = df.append(df_diff)
sns.pointplot(x = "samle number", y = "sensitivity", hue = "difference", data = df)
plot2([0.5, 0.8, 1.2, 2.0])

通过上图,我们发现,即使标准差达到0.10,当样本数量达到30以上时,T检验的效果就非常好了。

那么如果两个总体没有差异,T检验得出有差异的错误结论的情况是什么样子的呢?我们构建两个均值为1容量为1000的总体,每次随机抽取3个样本,在不同标准差及样本数目的情况下,考察T检验得出正确结论的频率。请看代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
def norm_t2(loc, scale1, scale2, n):
normSpace1 = stats.norm.rvs(loc = loc, scale = scale1, size = 1000)
normSpace2 = stats.norm.rvs(loc = loc, scale = scale2, size = 1000)
count = 0
for i in range(1000):
sample1 = np.random.choice(normSpace1, n)
sample2 = np.random.choice(normSpace2, n)
t, p = stats.ttest_ind(sample1, sample2)
if p &gt; 0.05:
count = count + 1
return count / 1000
def plot3(stds):
stds = [0.1, 0.5, 1]
df = pd.DataFrame()
for std in stds:
x = list(range(3, 60, 3))
y = [norm_t2(loc = 1, scale1 = std , scale2 = std, n = i) for i in x]
z = [std for i in x]
df_std = pd.DataFrame([x, y, z]).T
df_std.columns = ["samle number", "specificity", "std"]
if df.shape == (0, 0):
df = df_std
else:
df = df.append(df_std)
sns.pointplot(x = "samle number", y = "specificity", hue = "std", data = df)
plt.ylim(0, 1.2)
plot3([0.1, 0.5, 1])

通过上图,我们惊讶地发现,当两个总体没有差异时,不管标准差是多少和取样数目是多少,T检验都有非常好的表现,即得到差异不显著的的结论。

这里构建的总体都是正太总体,其实对于其他分布的总体,结论是一样的,这里就不展示了。

通过上面的探索,我们可以得到如下结论:

  1. 对于两个未知总体,通过T检验考察均值是否有显著差异,如果得出结论是差异不显著,那么进一步分析这些数据的标准差是否太大了,考虑是否增加抽样样本数做进一步分析。
  2. 对于两个未知总体,通过T检验考察均值是否有显著差异,如果得出结论是差异显著,那么请相信它吧!

如果差异不显著,我们通过增大样本数量,使得差异显著;如果差异显著,那就是一个好结果嘛!
哈哈,T检验是实验生物学家的利器啊!