python生信小练习（一） python生信小练习（一）

出自同哥的小练习，用于巩固基础知识：

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

用到的知识点
split
字符串的索引
输出格式为：

【python生信小练习（一）】NM_001011874
gcggcggcgggcgagcgggcgctggagtaggagctg.......

Answer:

for line in open(r'E:\Bioinformatics\Python\practice\chentong\notebook-master\data\test2.fa'): if line.startswith('>'): print(line.split(' ')[0]) else: print(line.strip('\n'))

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

join
strip
用到的知识点
输出格式为:

NM_001011874
gcggcggcgggc......TCCGCTG......GCGTTCACC......CGGGGTCCGGAG

Answer:

fasta = {} for line in open(r'E:\Bioinformatics\Python\practice\chentong\notebook-master\data\test2.fa'): if line.startswith('>'): key = line.split(' ')[0]#首行用空格做分隔符，取第一个元素作为基因名 fasta[key] = []#建立以基因名为key的字典，准备存储序列 else: fasta[key].append(line.strip())#若首行第一个字符不为'>'，将该行去除分隔符后，作为’值‘添加入字典for key, value in fasta.items(): print(key) print(''.join(value))

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

字符串切片操作
range
用到的知识点
输出格式为

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

Answer:

fasta = {} length = 80 for line in open(r'E:\Bioinformatics\Python\practice\chentong\notebook-master\data\test2.fa'): if line.startswith('>'): key = line.split(' ')[0]#首行用空格做分隔符，取第一个元素作为基因名 fasta[key] = []#建立以基因名为key的字典，准备存储序列 else: fasta[key].append(line.strip())#若首行第一个字符不为'>'，将该行去除分隔符后，作为’值‘添加入字典for key, value in fasta.items(): print(key) fasta_2 = ''.join(value) for i in range(0, len(fasta_2), length):#range: 1st: start; 2nd: end; 3rd: step length print(fasta_2[i : i + length])#slicing operation

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

sort
dict
aDict[key] = []
aDict[key].append(value)

fasta = {} for line in open(r'E:\Bioinformatics\Python\practice\chentong\notebook-master\data\test2.fa'): if line.startswith('>'): key = line.split(' ')[0]#首行用空格做分隔符，取第一个元素作为基因名 fasta[key] = []#建立以基因名为key的字典，准备存储序列 else: fasta[key].append(line.strip())#若首行第一个字符不为'>'，将该行去除分隔符后，作为’值‘添加入字典keyS = sorted(fasta.keys())#sorted函数按key值对字典排序 for i in keyS: print(i) print(''.join(fasta[i]))

5.1 提取给定名字的序列
用到的知识点

print >>fh, or fh.write()
取模运算，4 % 2 == 0
写程序 grepFasta.py, 提取fasta.name中名字对应的test2.fa的序列，并输出到屏幕。

Answer:

fasta = {} for line in open(r'E:\Bioinformatics\Python\practice\chentong\notebook-master\data\test2.fa'): if line.startswith('>'): key = line.split(' ')[0]#首行用空格做分隔符，取第一个元素作为基因名 fasta[key] = []#建立以基因名为key的字典，准备存储序列 else: fasta[key].append(line.strip())#读取fasta.name文件，并将其存储在name变量中 name = [] for line in open(r'E:\Bioinformatics\Python\practice\chentong\notebook-master\data\fasta.name'): line = '>' + line name.append(line.strip())#输出fasta[name]序列至屏幕 for i in name: print(i) print(''.join(fasta[i]))

5.2 写程序 grepFastq.py, 提取fastq.name中名字对应的test1.fq的序列，并输出到文件。
Answer:

fastq = {} i = 1 for line in open(r'E:\Bioinformatics\Python\practice\chentong\notebook-master\data\test1.fq'): if i % 4 == 1:#每条reads信息的第一行为reads名称 seqID = line.strip('\n')[1:]#取@之后的部分作为reads名 fastq[seqID] = [] elif i % 4 == 2: fastq[seqID] = line.strip('\n') i += 1#写入文件，出了一些问题，需要debug... FUCK!!!!!!!!! with open(r'E:\Bioinformatics\Python\practice\chentong\notebook-master\data\fastq.required') as f: for line in open(r'E:\Bioinformatics\Python\practice\chentong\notebook-master\data\fastq.name'): name = line.strip() f.write(str(name)) f.write(str(''.join(fastq[name])))

写程序 screenResult.py, 筛选test.expr中foldChange大于2的基因并且padj小于0.05的基，可以输出整行或只输出基因名字

逻辑与操作符 and
文件中读取的内容都为字符串，需要用int转换为整数，float转换为浮点数
Answer:

#Linux下 awk 方法，很简单： awk '{if($10 > 2 && $13 < 0.05) print $1}' ./test.expr | less#仅显示基因名称

文章图片
awk方法，仅显示基因名称

#R 语言实现： > f <- read.delim("clipboard", header = T, stringsAsFactors = F) > f[which(f$foldChange > 2 & f$padj < 0.05), 1] [1] "Novel00011" "Novel00043" "Novel00047" "Novel00077" "Novel00079" [6] "Novel00080" "Novel00084" "Novel00085" "Novel00086" "Novel00087" [11] "Novel00090" "Novel00124" "Novel00148" "Novel00156" "Novel00162" [16] "Novel00166"

#python实现： file = r'E:\Bioinformatics\Python\practice\chentong\notebook-master\data\test.expr' with open(file) as f: all = f.readlines()#读取文件所有行 lines = all[1:]#去除第一行 for item in lines: info = item.split('\t')#数据之间以换行符分割 if float(info[9].strip()) > 2.0 and float(info[12].strip()) < 0.05:#通过powershell以及debug看到， print(info[0])#每行后面的数字除了换行符还有空格，所以，这一步在转变数据类型之前进行空格的去除