python 处理生信数据基本操作

最新推荐文章于 2024-08-19 21:48:46 发布

qq_39306047

最新推荐文章于 2024-08-19 21:48:46 发布

阅读量1.6k

点赞数 1

分类专栏： python 文章标签：生物信息python

本文链接：https://blog.csdn.net/qq_39306047/article/details/90712707

版权

python 专栏收录该内容

58 篇文章 1 订阅

订阅专栏

作业1，2[给定FASTA格式的文件(test1.fa 和 test2.fa)，写一个程序 cat.py 读入文件，并输出到屏幕]代码：
for line in open(“c:/users/administer/desktop/test1.fa”):
print line.strip()

for line in open(“c:/users/administer/desktop/test1.fa”):
print (line.strip())

fjdjfkdjfkdsjfkdf

sdfjisdjfdj

sdfjkdjfkdjfkjkdjfk

sdkfjkdjfkjdskfjk777
&&kdjfjdfkdjfk
dkfjkdjfk
dfjidjfkdfj

c:/users/administer/desktop/test1.fa

for line in open(“c:/users/administer/desktop/test1.fa”):
print (line)

fjdjfkdjfkdsjfkdf

sdfjisdjfdj

sdfjkdjfkdjfkjkdjfk

sdkfjkdjfkjdskfjk777

&&kdjfjdfkdjfk

dkfjkdjfk

dfjidjfkdfj

c:/users/administer/desktop/test1.fa

3写程序 splitName.py, 读入test2.fa, 并取原始序列名字第一个空格前的名字为处理后的序列名字，输出到屏幕

split

字符串的索引

用到的知识点

输出格式为：

NM_001011874
gcggcggcgggcgagcgggcgctggagtaggagctg…

#作业3
filename = “c:/users/administer/desktop/test1.fa”
for line in open(filename):
if line[0] == ‘>’:
lineL = line.split(’ ')
print (lineL[0])
else:
print (line)

4写程序 formatFasta.py, 读入test2.fa，把每条FASTA序列连成一行然后输出

join

strip

用到的知识点

输出格式为:

NM_001011874
gcggcggcgggc…TCCGCTG…GCGTTCACC…CGGGGTCCGGAG

#作业4
filename = “c:/users/administer/desktop/test1.fa”
aList = []
for line in open(filename):
if line[0]==’>’:
lineL=line.split(’ ‘)
if aList:
print (’’.join(aList))
print (‘youyong妈’)
aList=[]
name=lineL[0]
print (name)
else:
aList.append(line.strip())####条件不成了首先执行这一步

#不要忘了最后一个序列
print (’’.join(aList))

5写程序 formatFasta-2.py, 读入test2.fa，把每条FASTA序列分割成20个字母一行的序列

字符串切片操作

range

用到的知识点

输出格式为

NM_001011874
gcggcggcgc.(60个字母).TCCGCTGACG #(每行80个字母)
acgtgctacg.(60个字母).GCGTTCACCC
ACGTACGATG(最后一行可不足80个字母)

#作业5
filename=“c:/users/administer/desktop/test1.fa”
length = 20
aList = []
for line in open(filename):
if line[0] == ‘>’:
lineL = line.split(’ ')
if aList:
seq = ‘’.join(aList)
len_seq = len(seq)
for i in range(0,len_seq,length):
print (seq[i:i+length])
aList = []
name = lineL[0]
print (name)
else:
aList.append(line.strip())
#不要忘了最后一个序列
seq = ‘’.join(aList)
len_seq = len(seq)
for i in range(0,len_seq,length):
print (seq[i:i+length])

6写程序 sortFasta.py, 读入test2.fa, 并取原始序列名字第一个空格前的名字为处理后的序列名字，排序后输出

sort

dict

aDict[key] = []

aDict[key].append(value)

用到的知识点

提取给定名字的序列

用到的知识点

print >>fh, or fh.write()

取模运算，4 % 2 == 0

#作业6
aDict = {}

filename=“c:/users/administer/desktop/test1.fa”

for line in open(filename):
if line[0] == ‘>’:
key = line.split()[0]
aDict[key] = []

else:
    aDict[key].append(line.strip())#文本以>号开头才成立

#--------END reading--------------------
keyL = aDict.keys()

keyL.sort()###没这个属性所以用不了

###This is because in Python 3 getting the keys() of a Dict returns a dict_keys object instead of a list. Change this:
targets.keys()
to this:
list(targets.keys())###
修改：
keyL=sorted(list(keyL))

for key in keyL:
print ("%s\n%s" % (key, ‘’.join(aDict[key])))

7.1写程序 grepFasta.py, 提取nameFile.txt中名字对应的test2.fa的序列，并输出到屏幕。

#作业7.1
aDict = {}

seqFile==“NM_01gggggggjjjjtest1.fa”
nameFile=“c:/users/administer/desktop/nameFile.txt”

for line in open(seqFile):
if line[0] == ‘>’:
key = line.split()[0][1:] #注意去掉开头的’>’，只取序列的名字
aDict[key] = [] #字典的值是一个列表
else:
aDict[key].append(line.strip())
#--------END reading--------------------
for line in open(nameFile):
name = line.strip()
print (">%s\n%s" % (name, ‘’.join(aDict[name])))

7.2写程序 grepFastq.py, 提取nameFile.txt中名字对应的test2.fq的序列，并输出到文件。

#作业7.2
[

0%4
0

1%4
1

2%4
2

3%4
3

4%4
0

5%4
1

count=0
for i in [1,2,3,4,5,6]:
count +=1
print (count)

1
2
3
4
5
6
count=0
for line in open(seqFile):
count +=1
if count%4==1:
name=line.strip()[1:]
print (count)
print (name)

aDict={}
count=0
for line in open(seqFile):
count +=1
if count%4==1:
name=line.strip()[1:]

	aDict[name]=["ggg"]
	aDict

{‘NM_0112835 gene=Rp1 CDS=128-6412’: [‘ggg’]}
{‘NM_0112835 gene=Rp1 CDS=128-6412’: [‘ggg’], ‘NM_0112888 gene=Rp1 CDS=128-6412’: [‘ggg’]}
{‘NM_0112835 gene=Rp1 CDS=128-6412’: [‘ggg’], ‘NM_0112888 gene=Rp1 CDS=128-6412’: [‘ggg’], ‘NM_01128777 gene=Rp1 CDS=128-6412’: [‘ggg’]}

nameFile.txt
NM_0112835 gene=Rp1 CDS=128-6412
errrr

]

seqFile = “c:/users/administer/desktop/test2.fq”
nameFile = “c:/users/administer/desktop/nameFile.txt”

aDict = {}
count = 0
for line in open(seqFile):
count += 1
if count % 4 == 1:
name = line.strip()[1:] #去掉开头的@
aDict[name] = [line] #存储名字行
else:
aDict[name].append(line)
#-------END reading----------这行知识注释，可有可无------
for line in open(nameFile):
name = line.strip()
print (’’.join(aDict[name]))

8写程序 screenResult.py, 筛选test.expr中foldChange大于2的基因并且padj小于0.05的基因

逻辑与操作符 and

文件中读取的内容都为字符串，需要用int转换为整数，float转换为浮点数

用到的知识点%%%%%%%可以用Linux实现一下这个功能
[
test.expr
foldChange padj
gene1 1 0.03
gene2 3 0.04
gene3 5 0.046
gene4 7 0.049
gene5 2 0.99

head=1
for line in open(file):
if head:
head -=1
print(head)
continue
lineL=line.strip()#strip()是处理字符的需要修改为split()处理单词字段的
foldChange = float(lineL[6])
adjp = float(lineL[8])
if foldChange>2 and adjp<0.05:
print (lineL[0])

0
g
g
g

head=1
for line in open(file):
if head:
head -=1
print(head)
continue
lineL=line.split()
foldChange = float(lineL[1])
adjp = float(lineL[2])
if foldChange>2 and adjp<0.05:
print (lineL[0])

0
gene2
gene3
gene4

]

#作业8
file = “c:/users/administer/desktop/test.expr”

head = 1 #This is used to indicate there is one head line needing to skip
for line in open(file):
if head:
head -= 1
continue
#-----Begin data lines------
lineL = line.split()

foldChange = float(lineL[7])
adjP = float(lineL[9])
if foldChange > 2 and adjP < 0.05:
print lineL[0]

9写程序 transferMultipleColumToMatrix.py 将文件(multipleColExpr.txt)中基因在多个组织中的表达数据转换为矩阵形式

aDict[‘key’] = {}

aDict[‘key’][‘key2’] = value

if key not in aDict

aDict = {‘ENSG00000000003’: {“A-431”: 21.3, “A-549”, 32.5,…},”ENSG00000000003”:{},}

用到的知识点

输入格式(只需要前3列就可以)

Gene Sample Value Unit Abundance
ENSG00000000003 A-431 21.3 FPKM Medium
ENSG00000000003 A-549 32.5 FPKM Medium
ENSG00000000003 AN3-CA 38.2 FPKM Medium
ENSG00000000003 BEWO 31.4 FPKM Medium
ENSG00000000003 CACO-2 63.9 FPKM High
ENSG00000000005 A-431 0.0 FPKM Not detected
ENSG00000000005 A-549 0.0 FPKM Not detected
ENSG00000000005 AN3-CA 0.0 FPKM Not detected
ENSG00000000005 BEWO 0.0 FPKM Not detected
ENSG00000000005 CACO-2 0.0 FPKM Not detected
输出格式

Name A-431 A-549 AN3-CA BEWO CACO-2
ENSG00000000460 25.2 14.2 10.6 24.4 14.2
ENSG00000000938 0.0 0.0 0.0 0.0 0.0
ENSG00000001084 19.1 155.1 24.4 12.6 23.5
ENSG00000000457 2.8 3.4 3.8 5.8 2.9

【>>> head = 1

tissueSet = set()
aDict = {}
for line in open(file):
if head:
head -=1
continue
lineL=line.split()
gene=lineL[0]

tissue=lineL[1]
expr=lineL[2]
if gene not in aDict:
	aDict[gene]={}
	aDict
assert tissue not in aDict[gene],"Duplicate tissue"
print (tissue)

{‘ENSG00000000003’: {}}
A-431
A-549
AN3-CA
BEWO
CACO-2
{‘ENSG00000000003’: {}, ‘ENSG00000000005’: {}}
A-431
A-549
AN3-CA
BEWO
CACO-2

for line in open(file):
if head:
head -=1
continue
lineL=line.split()
gene=lineL[0]

tissue=lineL[1]
expr=lineL[2]
if gene not in aDict:
	aDict[gene]={}
	
assert tissue not in aDict[gene],"Duplicate tissue"
aDict[gene][tissue]=expr
tissueSet.add(tissue)
aDict

{‘ENSG00000000003’: {‘A-431’: ‘21.3’}}
{‘ENSG00000000003’: {‘A-431’: ‘21.3’, ‘A-549’: ‘32.5’}}
{‘ENSG00000000003’: {‘A-431’: ‘21.3’, ‘A-549’: ‘32.5’, ‘AN3-CA’: ‘38.2’}}
{‘ENSG00000000003’: {‘A-431’: ‘21.3’, ‘A-549’: ‘32.5’, ‘AN3-CA’: ‘38.2’, ‘BEWO’: ‘31.4’}}
{‘ENSG00000000003’: {‘A-431’: ‘21.3’, ‘A-549’: ‘32.5’, ‘AN3-CA’: ‘38.2’, ‘BEWO’: ‘31.4’, ‘CACO-2’: ‘63.9’}}
{‘ENSG00000000003’: {‘A-431’: ‘21.3’, ‘A-549’: ‘32.5’, ‘AN3-CA’: ‘38.2’, ‘BEWO’: ‘31.4’, ‘CACO-2’: ‘63.9’}, ‘ENSG00000000005’: {‘A-431’: ‘0.0’}}
{‘ENSG00000000003’: {‘A-431’: ‘21.3’, ‘A-549’: ‘32.5’, ‘AN3-CA’: ‘38.2’, ‘BEWO’: ‘31.4’, ‘CACO-2’: ‘63.9’}, ‘ENSG00000000005’: {‘A-431’: ‘0.0’, ‘A-549’: ‘0.0’}}
{‘ENSG00000000003’: {‘A-431’: ‘21.3’, ‘A-549’: ‘32.5’, ‘AN3-CA’: ‘38.2’, ‘BEWO’: ‘31.4’, ‘CACO-2’: ‘63.9’}, ‘ENSG00000000005’: {‘A-431’: ‘0.0’, ‘A-549’: ‘0.0’, ‘AN3-CA’: ‘0.0’}}
{‘ENSG00000000003’: {‘A-431’: ‘21.3’, ‘A-549’: ‘32.5’, ‘AN3-CA’: ‘38.2’, ‘BEWO’: ‘31.4’, ‘CACO-2’: ‘63.9’}, ‘ENSG00000000005’: {‘A-431’: ‘0.0’, ‘A-549’: ‘0.0’, ‘AN3-CA’: ‘0.0’, ‘BEWO’: ‘0.0’}}
{‘ENSG00000000003’: {‘A-431’: ‘21.3’, ‘A-549’: ‘32.5’, ‘AN3-CA’: ‘38.2’, ‘BEWO’: ‘31.4’, ‘CACO-2’: ‘63.9’}, ‘ENSG00000000005’: {‘A-431’: ‘0.0’, ‘A-549’: ‘0.0’, ‘AN3-CA’: ‘0.0’, ‘BEWO’: ‘0.0’, ‘CACO-2’: ‘0.0’}}
】
#作业9

file = “c:/users/administer/desktop/multipleColExpr.txt”
head = 1
tissueSet = set()
aDict = {}
for line in open(file):
if head: #skip header line
head -= 1
continue
#-------------------------------
lineL = line.split()
gene = lineL[0]
tissue = lineL[1]
expr = lineL[2]
if gene not in aDict:
aDict[gene] = {}
assert tissue not in aDict[gene], “Duplicate tissues”
aDict[gene][tissue] = expr
tissueSet.add(tissue)
#—END reading----------------------------
tissueL = list(tissueSet)
tissueL.sort()
print (“Gene\t%s” % ‘\t’.join(tissueL))

for gene, tissueD in aDict.items():
exprL = [gene]
for tissue in tissueL:
exprL.append(tissueD[tissue])
print (’\t’.join(exprL))
Overwriting transferMultipleColumToMatrix.py
%run script/transferMultipleColumToMatrix
Gene A-431 A-549 AN3-CA BEWO CACO-2
ENSG00000000460 25.2 14.2 10.6 24.4 14.2
ENSG00000000938 0.0 0.0 0.0 0.0 0.0
ENSG00000000457 2.8 3.4 3.8 5.8 2.9
ENSG00000000419 73.8 38.6 33.9 53.7 155.5
ENSG00000000003 21.3 32.5 38.2 31.4 63.9
ENSG00000000005 0.0 0.0 0.0 0.0 0.0

10写程序 reverseComplementary.py计算序列 ACGTACGTACGTCACGTCAGCTAGAC的反向互补序列

reverse

list(seq)

用到的知识点

#作业10
str1 = “ACGTACGTACGTCACGTCAGCTAGAC”

ATCG_dict = {‘A’:‘T’,‘T’:‘A’,‘C’:‘G’,‘G’:‘C’, ‘a’:‘t’,‘t’:‘a’,‘c’:‘g’,‘g’:‘c’}

tmpL = []
for i in str1:
tmpL.append(ATCG_dict[i])
tmpL.reverse()
print (’’.join(tmpL))

11写程序 collapsemiRNAreads.py转换smRNA-Seq的测序数据

输入文件格式(mir.collapse, tab-分割的两列文件，第一列为序列，第二列为序列被测到的次数)

ID_REF VALUE
ACTGCCCTAAGTGCTCCTTCTGGC 2
ATAAGGTGCATCTAGTGCAGATA 25
TGAGGTAGTAGTTTGTGCTGTTT 100
TCCTACGAGTTGCATGGATTC 4
输出文件格式 (mir.collapse.fa, 名字的前3个字母为样品的特异标示，中间的数字表示第几条序列，是序列名字的唯一标示，
第三部分是x加每个reads被测到的次数。三部分用下划线连起来作为fasta序列的名字。)

ESB_1_x2
ACTGCCCTAAGTGCTCCTTCTGGC
ESB_2_x25
ATAAGGTGCATCTAGTGCAGATA
ESB_3_x100
TGAGGTAGTAGTTTGTGCTGTTT
ESB_4_x4
TCCTACGAGTTGCATGGATTC

#作业11
smRNA_file = “c:/users/administer/desktop/mir.collapse”

head = 1
sample = ‘ESB’
lineno = 0
for line in open(smRNA_file):
if head:
head -= 1
continue
#----end skip header line—
seq, value = line.split()
lineno += 1
print (">%s_%d_x%s\n%s" % (sample, lineno, value, seq))

12简化的短序列匹配程序 (map.py) 把short.fa中的序列比对到ref.fa, 输出短序列匹配到ref.fa文件中哪些序列的哪些位置

find

用到的知识点

输出格式 (输出格式为bed格式，第一列为匹配到的染色体，第二列和第三列为匹配到染色体序列的起始终止位置
（位置标记以0为起始，代表第一个位置；终止位置不包含在内，第一个例子中所示序列的位置是(199,208]
(前闭后开，实际是chr1染色体第199-206的序列，0起始). 第4列为短序列自身的序列.)。

附加要求：可以只匹配到给定的模板链，也可以考虑匹配到模板链的互补链。这时第5列可以为短序列的名字，
第六列为链的信息，匹配到模板链为’+’，匹配到互补链为’-‘。注意匹配到互补链时起始位置也是从模板链的5’端算起的。

chr1 199 208 TGGCGTTCA
chr1 207 216 ACCCCGCTG
chr2 63 70 AAATTGC
chr3 0 7 AATAAAT

#作业12
ref = “data/ref.fa”
short = “data/short.fa”

refD = {}
for line in open(ref):
if line[0] == ‘>’:
name = line[1:-1]
assert name not in refD
refD[name] = []####{chr1:[],…}
else:
refD[name].append(line.strip())
#---------END reading----------------
for key, valueL in refD.items():
refD[key] = ‘’.join(valueL)###chr1:kkkkkklianyiqi
#------------------------------------
for line in open(short):
if line[0] == ‘>’:
pass
else:
seq = line.strip()
len_seq = len(seq)
for key, value in refD.items():
pos = value.find(seq)
while pos != -1:
print “%s\t%d\t%d\t%s” % (key, pos, pos+len_seq, seq)
current = pos + 1
newpos = value[current:].find(seq)
if newpos == -1:
break
pos = current + newpos

Overwriting map.py
%run script/map
chr1 199 208 TGGCGTTCA
chr1 207 216 ACCCCGCTG
chr2 63 70 AAATTGC
chr3 0 7 AATAAAT
chr3 66 74 CCACACAA
chr3 206 214 ATATCACT
chr2 163 171 ATATCACT
chr3 417 423 AGAGGA
chr2 374 380 AGAGGA
chrX 137 143 AGAGGA