病毒基因分类

输入多个国家的新型冠状病毒的基因组数据，通过子串（kmer）词频分析，挖掘国家间的特征子串，用于病毒的分类。

1.数据预处理

（1）字符串分割与子串统计

def kmer(str1,N):"""str1：要分割的字符串N:子串的长度return :返回一个统计好的字典"""import copy#存放子串str_dict = {}#子串长度#N = 5#深度复制str2 = copy.deepcopy(str1)#去掉首尾符号，并转换为liststr2 = list(str2.strip())#str2长度小于子串长度时子串切割完毕，退出循环while len(str2) >= N:str_tam = str2[0]for str3 in str2[1:N]:str_tam = str_tam + str3#若str_tam已存在str_dict中则加1，否则添加进str_dict中if str_tam in str_dict.keys():str_dict[str_tam] = str_dict[str_tam] + 1else:str_dict[str_tam] = 1#删除第一字母del str2[0]return str_dict

str1 = 'ABCDEFGHABCDEF'
str_dict = kmer(str1,5)

str_dict

{'ABCDE': 2,'BCDEF': 2,'CDEFG': 1,'DEFGH': 1,'EFGHA': 1,'FGHAB': 1,'GHABC': 1,'HABCD': 1}

（2）基因文件处理

def processing_file(filepath):"""基因文件处理"""import redata = []tam_list = []for line in open(filepath):#判断行是否为空if line.strip() is not '':#匹配，判断是字段信息还是基因序列if re.match(">", line):#初始化临时列表tam_list = []#去掉首尾字符，并分割字符串line_list = line.strip().split('/')#去掉">"得到病毒名称name = line_list[0].strip(">").strip()#国家country = line_list[1]#基因编号ids = line_list[2]#分割最后一个字符串tam = line_list[3].split("|")#时间date = tam[0]#tl = tam[1]#添加到临时列表tam_list.append(name)tam_list.append(country)tam_list.append(ids)tam_list.append(date)tam_list.append(tl)else:#基因序列line = line.strip()#添加到临时列表tam_list.append(line)#把临时列表添加到data列表data.append(tam_list)return data

filepath = 'test.fa'

data = processing_file(filepath)

len(data)

20

2.存入MySQL数据库

https://www.runoob.com/python3/python3-mysql.html

import pymysql
#连接数据库
def connectdb():# 打开数据库连接db = pymysql.connect(host='127.0.0.1', port=3306, db='bingdubase', user='root', passwd='*****', charset='utf8')return db#创建表
def clear_table(db,table):# 使用cursor()方法获取操作游标 cursor = db.cursor()# 使用 execute() 方法执行 SQL，如果表存在则删除cursor.execute("DROP TABLE IF EXISTS "+table)# 使用预处理语句创建表sql = """CREATE TABLE """+table+"""(name VARCHAR(255) NOT NULL,country VARCHAR(255) NOT NULL,ids VARCHAR(255) NOT NULL,date VARCHAR(255) NOT NULL,tl VARCHAR(255) NOT NULL,seq LONGTEXT NOT NULL,INCOME FLOAT )"""# 执行sql语句cursor.execute(sql)# 执行sql语句db.commit()#插入数据库
def insert(db,data):# 使用cursor()方法获取操作游标 cursor = db.cursor()# SQL 插入语句sql = "INSERT INTO gene_info(name,country,ids,date,tl,seq)\VALUES('%s','%s','%s','%s','%s','%s')"%(data[0],data[1],data[2],data[3],data[4],data[5])# 执行sql语句cursor.execute(sql)# 执行sql语句db.commit()#查询
def select(db,table):# 使用cursor()方法获取操作游标 cursor = db.cursor()#查询sql语句sql = "SELECT * FROM "+table# 执行SQL语句cursor.execute(sql)# 获取所有记录列表results = cursor.fetchall()return results
#查询
def select_re_ids(db,table):# 使用cursor()方法获取操作游标 cursor = db.cursor()#查询sql语句sql = "SELECT ids FROM "+table# 执行SQL语句cursor.execute(sql)# 获取所有记录列表results = cursor.fetchall()return results

def save_data_to_mysql():"""存数据到mysql数据库"""filepath = "test.fa"#处理基因数据data = processing_file(filepath)#存到数据库#连接数据库db = connectdb()#建表clear_table(db,"gene_info")for line in data:try:insert(db,line)except:passdb.close()

#save_data_to_mysql()

3.kmer预处理并存到Hbase数据库

Python读取Mysql数据，并对基因序列做kmer预处理，参考下列结构存储到MongoDB数据库或者Hbase数据库。

{Seqid:” ”, Country:” ”, kmer1:number, kmer2:number,-------- }

https://blog.csdn.net/y472360651/article/details/79059457

import happybase
#连接hbase
def connHbase():conn = happybase.Connection(host="127.0.0.1",port=9090)return conn

#关闭连接
#conn.close()

create_table(name,families)：创建表，无返回值

name：表名

families：列族

def createTable(tablename):"""建表"""#tablename='gene_info'conn = connHbase()families = {"baseinfo":dict(),"kmerinfo":dict()}conn.create_table(tablename,families)      # 如果连接时，有传递表前缀参数时，真实表名将会是："{}_{}".format(table_prefix,name)conn.close()

delete_table(name,disable=False)：删除表，无返回值

name：表名

disable：是否先禁用表

def deleteTable(tablename):conn = connHbase()#先禁用表conn.disable_table(tablename)conn.delete_table(tablename,disable=False)

表(happybase.Table)

happybase.Table(name,connection)：获取表实例

name：表名

connection：连接

put(row, data, timestamp=None, wal=True)：插入数据，无返回值

row: 行

data: 数据，dict类型，{列:值}构成，列与值皆为str类型

timestamp：时间戳，默认None，即写入当前时间戳

wal：是否写入wal，默认为True

def hbase_main():tablename = 'gene_info'#连接hbaseconn = connHbase()#获取所有表able_name_list = conn.tables()conn.close()#判断表是否存在if b'gene_info' in able_name_list:#删除表deleteTable(tablename)#创建表createTable(tablename) #数据kmer处理并存入hbase#tablename = "gene_info"#子串长度N = 5db = connectdb()results = select(db,tablename)db.close()#连接hbaseconn = connHbase()if conn.is_table_enabled(tablename):passelse:#启用表conn.enable_table(tablename)#连接表table = happybase.Table(tablename,conn)info = table.families()for row in results:name = str(row[0])country = str(row[1])ids = str(row[2])date = str(row[3])tl = str(row[4])table.put(ids,{"baseinfo:name":name})table.put(ids,{"baseinfo:country":country})table.put(ids,{"baseinfo:date":date})table.put(ids,{"baseinfo:tl":tl})#kmer处理seq = row[5]data_list = kmer(seq,N)for key in data_list.keys():table.put(ids,{"kmerinfo:"+key:str(data_list[key])})conn.close()print("已存入Hbase数据库")

#hbase_main()

4.kmer进行降维处理

Python读取Hbase数据，并对kmer进行降维处理，找到类间频率差异最大的n个kmer子串.

https://www.runoob.com/numpy/numpy-matrix.html

降维：https://www.cnblogs.com/yanqiang/p/11781377.html

降维处理：归一化

将样本点的数值减去最小值，再除以样本点数值最大与最小的差

#归一化处理
def noramlization(data):data = np.array(data)minVals = data.min(0)maxVals = data.max(0)ranges = maxVals - minValsm = data.shape[0]normData = data - minValsnormData = normData/rangesreturn normData

tablename = "gene_info"
db = connectdb()
results = select_re_ids(db,tablename)
db.close()

#获取行号
rows = []
for ids in results:rows.append(ids[0])

columns = ["baseinfo:country","kmerinfo"]
conn = connHbase()
table = happybase.Table(tablename,conn)
scanner = table.rows(rows=rows,columns=columns,timestamp=None, include_timestamp=False)
conn.close()

len(scanner)

20

X = []  #存特征值
Y = []  #存标签：国家
X_len_min = 2000
for info in scanner:#获取kmer的值，并转换成listvalues_list = list(info[1].values())country = values_list[0]values = []for value in values_list[1:]:#bytes 转换为 intvalues.append(int(value))if X_len_min > len(values):X_len_min = len(values)X.append(values)Y.append(country)

X_len_min

1028

#数据统一长度
X_new = []
for xn in X:X_new.append(xn[:X_len_min])

len(X_new[14])

1028

import numpy as np
import matplotlib.pyplot as plt#归一化处理
def noramlization(data):minVals = data.min(0)maxVals = data.max(0)ranges = maxVals - minValsnormData = np.zeros(np.shape(data))m = data.shape[0]normData = data - np.tile(minVals, (m, 1))normData = normData/np.tile(ranges, (m, 1))return normData, ranges, minValsx = np.array(X_new)X, _, _ = noramlization(x)

5.分类

import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
import warnings
from sklearn import svm#svm导入
from sklearn.svm import SVC
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from sklearn.exceptions import ChangedBehaviorWarning## 设置属性防止中文乱码
mpl.rcParams['font.sans-serif'] = [u'SimHei']
mpl.rcParams['axes.unicode_minus'] = Falsewarnings.filterwarnings('ignore', category=ChangedBehaviorWarning)x, y = X,Y
y = pd.Categorical(y).codes #把文本数据进行编码，比如a b c编码为 0 1 2
## 数据分割
x_train, x_test, y_train, y_test = train_test_split(x, y, random_state=0, train_size=0.8)

使用SVM分类器进行模型构建

clf = SVC(C=1,kernel='rbf',gamma=0.1)## 模型训练
clf.fit(x_train, y_train)

SVC(C=1, break_ties=False, cache_size=200, class_weight=None, coef0=0.0,decision_function_shape='ovr', degree=3, gamma=0.1, kernel='rbf',max_iter=-1, probability=False, random_state=None, shrinking=True,tol=0.001, verbose=False)

计算模型的准确率/精度

print (clf.score(x_train, y_train)) 
print ('训练集准确率：', accuracy_score(y_train, clf.predict(x_train)))
print (clf.score(x_test, y_test))
print ('测试集准确率：', accuracy_score(y_test, clf.predict(x_test)))

1.0
训练集准确率： 1.0
0.5
测试集准确率： 0.5

计算决策函数的结构值以及预测值(decision_function计算的是样本x到各个分割平面的距离<也就是决策函数的值>)

print ('decision_function:\n', clf.decision_function(x_train))
print ('\npredict:\n', clf.predict(x_train))

decision_function:[ 0.9993134   0.99991139 -1.00048851  0.99933992 -1.00010121 -0.99985233-0.99985233  1.00014379 -0.99985621 -0.99971422  1.00028099  0.99932378-1.00060865 -1.00006369  0.99952865  0.99934374]predict:[1 1 0 1 0 0 0 1 0 0 1 1 0 0 1 1]