生信算法1 - DNA测序算法实践之序列操作
2023-12-13 14:49:28
生信序列基本操作算法
建议在Jupyter实践,python版本3.9
1. 定义DNA序列
seq = 'ACGT'
# 从序列指定索引的碱基
seq[1]
# 'C'
# 序列长度
len(seq)
# 4
2. 序列拼接
# 序列拼接 - 字符串
seq1 = 'AACC'
seq2 = 'GGTT'
print(seq1 + seq2)
# AACCGGTT
# 序列拼接 - 列表
seqs = ['A', 'C', 'G', 'T']
print(''.join(seqs))
3. 序列生成
# 随机产生ATCG碱基
import random
random.choice('ACGT')
# 生成DNA序列
import random
seq = ''
for _ in range(100):
seq += random.choice('ACGT')
print(seq)
# GTGTCGACACCCGAGTGTTCACAAGTGGTCTTTCGGCACTTTGCTCTGAGTGGCGCGCTTAGCTATAAGGGCAAGTATCCTAGTCATATCCCGAAGGCGT
# 生成DNA序列方法2
seq = ''.join([random.choice('ACGT') for _ in range(10)])
print(seq)
# ATGGATGTCT
4. 序列截取
seq = "ATGGATGTCT"
seq[1:3]
# 'AT'
seq[:3]
#'GAT'
seq[7:]
# 'GAC'
seq[-3:]
# 'GAC'
5. 序列比较
# 获取序列最长相同前缀
def longestCommonPrefix(s1, s2):
i = 0
while i < len(s1) and i < len(s2) and s1[i] == s2[i]:
i += 1
return s1[:i]
longestCommonPrefix('ACCATTG', 'ACCAAGTC')
# 'ACCA'
# 比较序列是否相同
def match(s1, s2):
if not len(s1) == len(s2):
return False
for i in range(0, len(s1)):
if not s1[i] == s2[i]:
return False
return True
match('ACCATTG', 'ACCATTG')
# False
# 序列字符串比较
'ACCATTG' == 'ACCATTG'
True
5. 获取互补序列
def Complement(seq):
complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
t = ''
for base in seq:
t = t + complement[base]
return t
Complement('ACCATTG')
# 'TGGTAAC'
6. 获取反向互补序列
def reverseComplement(seq):
complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
t = ''
for base in seq:
t = complement[base] + t
return t
reverseComplement('ACCATTG')
# 'CAATGGT'
7. 读取基因组fasta序列文件
def readGenome(filename):
genome = ''
with open(filename, 'r') as f:
for line in f:
# 忽略文件头信息
if not line[0] == '>':
genome += line.rstrip()
return genome
genome = readGenome('lambda_virus.fa')
len(genome)
# 48502
genome[:100]
# 'GGGCGGCGACCTCGCGGGTTTTCGCTATTTATGAAAATTTTCCGGTTTAAGGCGTTTCCGTTCTTCTTCGTCATAACTTAATGTTTTTATTTAAAATACC'
7. 计算基因组fasta文件ATCG碱基的数量
# 循环遍历方法
counts = {'A': 0, 'C': 0, 'G': 0, 'T': 0}
for base in genome:
counts[base] += 1
print(counts)
# {'A': 12334, 'C': 11362, 'G': 12820, 'T': 11986}
# 调用库
import collections
collections.Counter(genome)
# Counter({'G': 12820, 'A': 12334, 'T': 11986, 'C': 11362})
文章来源:https://blog.csdn.net/LittleComputerRobot/article/details/134969741
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。 如若内容造成侵权/违法违规/事实不符,请联系我的编程经验分享网邮箱:veading@qq.com进行投诉反馈,一经查实,立即删除!
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。 如若内容造成侵权/违法违规/事实不符,请联系我的编程经验分享网邮箱:veading@qq.com进行投诉反馈,一经查实,立即删除!