生物信息中的Python 01 | 从零开始处理基因序列
一、 序列数据的下载
在开始了解序列的处理流程时，我们先要知道序列下载网址。其中一个知名的网站就是NCBI （National Center for Biotechnology Information）美国国立生物技术信息中心。
1、通过如下的网站进入 NCBI ，可以看到它包含许多的子库，其中 Gene 就是我们一般下载基因序列的库，接下来，在后面的输入框输入 oct4 并点击 Search。
NCBI: https://www.ncbi.nlm.nih.gov/
2、可以看到该基因在不同物种和实验中所测得的相同基因序列，我们选择其中智人的POU5F1基因。
值得注意的是 POU5F1 是 Oct4 基因的别名，本质上指的一个基因
3、向下滚动，直到看到如下图所示的 FASTA 链接，点击进入。
4、在这个页面就可以看到通过测序技术所得到的DNA序列。
5、通过如下步骤我们可以得到该基因序列的 fasta 格式文件
6、你也可以按照上述步骤尝试获取[ Mus musculus ] 的 fasta 序列，我们后面的分析需要用到
二、 DNA序列基本处理
Python版本：Python 3.6
IDE：Pycharm (https://www.jetbrains.com/pycharm/) 下载 Pycharm 的免费社区版就足够我们学习使用
操作系统：Win7
0、在Pycharm里新建如下目录的项目
1、进入main.py文件，我们先把序列文件读取出来看看，到底是怎样的结果
with open('res/sequence1.fasta') as file:for line in file:print (line)
2、可以看到Fasta格式开始于一个标识符>，然后是一行描述，下面是序列，直到下一个>,表示下一条序列
这些字符串看起来和下载 Fasta 文件页面显示的差不多，但是这不是我们想要的结果
Fasta 格式详解
3、接下来我们把描述字段和序列分别提取并存储在字典中
fasta = {}
with open('res/sequence1.fasta') as file:sequence = ""for line in file:if line.startswith(">"):# 去除描述字段行中的\n和>name = line[1:].rstrip()fasta[name] = ''continue# 去除序列字段行中的\n，并将所有字符规范为大写字符fasta[name] += line.rstrip().upper()
print (fasta)
用函数把上面的代码装起来，方便后续调用
def get_fasta(fasta_path):fasta = {}with open(fasta_path) as file:sequence = ""for line in file:if line.startswith(">"):# 去除描述字段行中的\n和>name = line[1:].rstrip()fasta[name] = ''continue# 去除序列字段行中的\n，并将所有字符规范为大写字符fasta[name] += line.rstrip().upper()return fasta
4、拿到规范化的数据，我们现在来看看具有它具有的生物学意义，这里为了以后方便调用，使用函数的形式来实现
4.1 核苷酸计数，碱基偏好性：
这里的统计数值可以查看碱基偏好性。比如， 一定类型的小RNA会有特定的碱基偏好性，它的第一个碱基偏好U。可以用于评价数据质量。如果miRNA 第一碱基不是U偏好，说明数据或分析过程有问题。
# 核苷酸计数
def nt_count(seq):ntCounts = []for nt in ['A', 'C', 'G', 'T']:ntCounts.append(seq.count(nt))return ntCounts
4.2 GC含量：
（A+T）/（G+C）之比随DNA的种类不同而异。GC含量愈高，DNA的密度也愈高，同时热及碱不易使之变性，因此利用这一特性便可进行DNA的分离或测定。同时，物种的GC含量有着特异性，以此可以判断测序后的数据是否合格。
# CG 含量
from __future__ import division
def cg_content(seq):total = len(seq)gcCount = seq.count('G') + seq.count('C')gcContent = format(float(gcCount / total * 100), '.6f')return gcContent
4.3 DNA 翻译为 RNA：
# DNA 翻译为 RNA
def dna_trans_rna(seq):rnaSeq = re.sub('T', 'U', seq)# method2: rnaSeq = dnaSeq.replace('T', 'U')return rnaSeq
4.4 RNA 翻译为 蛋白质：
def rna_trans_protein(rnaSeq):codonTable = {'AUA':'I', 'AUC':'I', 'AUU':'I', 'AUG':'M','ACA':'T', 'ACC':'T', 'ACG':'T', 'ACU':'T','AAC':'N', 'AAU':'N', 'AAA':'K', 'AAG':'K','AGC':'S', 'AGU':'S', 'AGA':'R', 'AGG':'R','CUA':'L', 'CUC':'L', 'CUG':'L', 'CUU':'L','CCA':'P', 'CCC':'P', 'CCG':'P', 'CCU':'P','CAC':'H', 'CAU':'H', 'CAA':'Q', 'CAG':'Q','CGA':'R', 'CGC':'R', 'CGG':'R', 'CGU':'R','GUA':'V', 'GUC':'V', 'GUG':'V', 'GUU':'V','GCA':'A', 'GCC':'A', 'GCG':'A', 'GCU':'A','GAC':'D', 'GAU':'D', 'GAA':'E', 'GAG':'E','GGA':'G', 'GGC':'G', 'GGG':'G', 'GGU':'G','UCA':'S', 'UCC':'S', 'UCG':'S', 'UCU':'S','UUC':'F', 'UUU':'F', 'UUA':'L', 'UUG':'L','UAC':'Y', 'UAU':'Y', 'UAA':'', 'UAG':'','UGC':'C', 'UGU':'C', 'UGA':'', 'UGG':'W',}proteinSeq = ""for codonStart in range(0, len(rnaSeq), 3):codon = rnaSeq[codonStart:codonStart + 3]if codon in codonTable:proteinSeq += codonTable[codon]return proteinSeq
4.5 获取反向序列：
# 获取反向序列
def reverse_comple(type, seq):seq = seq[::-1]dnaTable = {"A":"T", "T":"A", "C":"G", "G":"C"}rnaTable = {"A": "T", "U": "A", "C": "G", "G": "C"}res = ""if type == "dna":for ele in seq:if ele in seq:if type == "dna":res += dnaTable[ele]else:res += rnaTable[ele]return res
4.6 最后我们来一个main来把上面的函数统统运行一遍
if __name__ == '__main__':oct4 = get_fasta('res/sequence1.fasta')for name, sequence in oct4.items():print ("name: ", name)print ("sequence: ", sequence)print ("nt_count: ", nt_count(sequence))print ("cg_content: ", cg_content(sequence))rna = dna_trans_rna(sequence)print ("rna: ", rna)protein = rna_trans_protein(rna)print ("protein: ", protein)print ("reverse_comple: ", reverse_comple("dna", sequence))
部分结果如下：
了解基因结构
生物信息中的Python 02 | 用biopython解析序列
生物信息中的Python 03 | 自动化操作NCBI
生物信息中的Python 04 | 批量下载基因与文献
生物信息中的Python 01 | 从零开始处理基因序列相关推荐
- 生物信息中的Python 02 | 用biopython解析序列
上一篇文章生物信息中的Python 01 | 从零开始处理基因序列自己造轮子实现了序列的基础操作,但是在Python的世界里,一项工作只要重复的次数多了,那么一定就会有大神来开发相应的包来解决,这个包 ...
- python中序列和列表区别细菌真菌病毒_生物信息中的Python 02 | 用biopython解析序列...
上一篇文章生物信息中的Python 01 | 从零开始处理基因序列自己造轮子实现了序列的基础操作,但是在Python的世界里,一项工作只要重复的次数多了,那么一定就会有大神来开发相应的包来解决,这个包 ...
- 生物信息中的Python 05 | 从 Genbank 文件中提取 CDS 等其他特征序列
1 介绍 在基因结构分析或其他生物功能分析中会时常用到 CDS 序列,以及其他诸如 mRNA 序列,misc RNA序列等具有生物意义的序列片段.而NCBI 的基因库中已经包含有这些的信息,但是只有一 ...
- 利用python处理dna序列_利用Python编程提取基因组基因序列
生物技术. DOI: 10.16660/j.cnki.1674-098X.2019.11.141 利用Python编程提取基因组基因序列 ① 庞雪原 张婷婷 (东北农业大学生命科学学院 黑龙江哈尔滨 ...
- Blender中的Python脚本介绍学习教程
Blender中的Python脚本介绍学习教程 MP4 |视频:h264,1280×720 |音频:AAC,48000 Hz 语言:英语+中英文字幕(根据原英文字幕机译更准确)|大小解压后:1.63 ...
- python写一个游戏多少代码-使用50行Python代码从零开始实现一个AI平衡小游戏
集智导读: 本文会为大家展示机器学习专家 Mike Shi 如何用 50 行 Python 代码创建一个 AI,使用增强学习技术,玩耍一个保持杆子平衡的小游戏.所用环境为标准的 OpenAI Gym, ...
- c中调用python解释器
c中调用python解释器 非常的简单, Code #include <iostream> #include <string> #include "Python.h& ...
- Windows中使用Python和C/C++联合开发应用程序起步
Windows中使用Python和C/C++联合开发应用程序起步(Python/C API应用) 2009-10-01 23:27:38| 分类: 学习 |举报 |字号订阅 C/C++几乎无所不能, ...
- java执行python脚本_使用Runtime.getRuntime().exec()在java中调用python脚本
举例有一个Python脚本叫test.py,现在想要在Java里调用这个脚本.假定这个test.py里面使用了拓展的包,使得pythoninterpreter之类内嵌的编译器无法使用,那么只能采用ja ...
最新文章
- 1480. 一维数组的动态和
- java读取配置文件的几种方法
- Windows下Python环境搭建
- SVN提示：由于目标机器积极拒绝，无法连接 的解决方法
- mac+修改+ssh文件夹权限_linux ,mac连接, git pull error, chmod修改文件的权限/chown修改文件和目录的所有者...
- 遥感RS植被指数大全整理
- mysql的表创建_MYSQL表创建
- oracle 批量插入 mysql 区别,Oracle 和 mysql 的批量操作Sql语句 的区别
- React.js 小书 Lesson24 - PropTypes 和组件参数验证
- Keil MDK下载程序时的相关设置
- mybatis加载属性
- logrotate日志轮转配置文档
- Delphi IDE下载全地址
- vs2015 安装破解版Visual Assist X
- web功能测试方法大全——完整！全面！（纯干货，建议收藏哦~）
- 阿里云对象存储OSS
- 主码流和辅码流的区别
- 可能是最好的正则表达式的教程笔记了吧...
- 云队友丨快手估值越高，越像抖音
- ElementUI 整体页面布局