python生信小练习（一）

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

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

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

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'))

2. 写程序 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))

3. 写程序 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

4. 写程序 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])))

6. 写程序 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])    #每行后面的数字除了换行符还有空格，所以，这一步在转变数据类型之前进行空格的去除

python result