出自同哥的小練習(xí)，用于鞏固基礎(chǔ)知識：

1. 寫程序 splitName.py, 讀入test2.fa, 并取原始序列名字第一個(gè)空格前的名字為處理后的序列名字，輸出到屏幕

• 用到的知識點(diǎn)
• 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
• 用到的知識點(diǎn)
• 輸出格式為:

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]    #首行用空格做分隔符，取第一個(gè)元素作為基因名
        fasta[key] = []             #建立以基因名為key的字典，準(zhǔn)備存儲序列
    else:
        fasta[key].append(line.strip())   #若首行第一個(gè)字符不為'>'，將該行去除分隔符后，作為’值‘添加入字典

for key, value in fasta.items(): 
    print(key)
    print(''.join(value))

3. 寫程序 formatFasta-2.py, 讀入test2.fa，把每條FASTA序列分割成80個(gè)字母一行的序列

• 字符串切片操作
• range
• 用到的知識點(diǎn)
• 輸出格式為

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

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]    #首行用空格做分隔符，取第一個(gè)元素作為基因名
        fasta[key] = []             #建立以基因名為key的字典，準(zhǔn)備存儲序列
    else:
        fasta[key].append(line.strip())   #若首行第一個(gè)字符不為'>'，將該行去除分隔符后，作為’值‘添加入字典

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, 并取原始序列名字第一個(gè)空格前的名字為處理后的序列名字，排序后輸出
   用到的知識點(diǎn)

• 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]    #首行用空格做分隔符，取第一個(gè)元素作為基因名
        fasta[key] = []             #建立以基因名為key的字典，準(zhǔn)備存儲序列
    else:
        fasta[key].append(line.strip())   #若首行第一個(gè)字符不為'>'，將該行去除分隔符后，作為’值‘添加入字典

keyS = sorted(fasta.keys())     #sorted函數(shù)按key值對字典排序
for i in keyS:
    print(i)
    print(''.join(fasta[i]))

5.1 提取給定名字的序列
用到的知識點(diǎn)

• print >>fh, or fh.write()
• 取模運(yùn)算，4 % 2 == 0
• 寫程序 grepFasta.py, 提取fasta.name中名字對應(yīng)的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]    #首行用空格做分隔符，取第一個(gè)元素作為基因名
        fasta[key] = []             #建立以基因名為key的字典，準(zhǔn)備存儲序列
    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中名字對應(yīng)的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
• 文件中讀取的內(nèi)容都為字符串，需要用int轉(zhuǎn)換為整數(shù)，float轉(zhuǎn)換為浮點(diǎn)數(shù)
  Answer:

#Linux下 awk 方法，很簡單：
awk '{if($10 > 2 && $13 < 0.05) print $1}' ./test.expr | less   #僅顯示基因名稱

awk方法，僅顯示基因名稱

#R 語言實(shí)現(xiàn)：
> 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實(shí)現(xiàn)：
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')   #數(shù)據(jù)之間以換行符分割
        if float(info[9].strip()) > 2.0 and float(info[12].strip()) < 0.05:    #通過powershell以及debug看到，
            print(info[0])    #每行后面的數(shù)字除了換行符還有空格，所以，這一步在轉(zhuǎn)變數(shù)據(jù)類型之前進(jìn)行空格的去除

python result