上一讲中我们初步认识了python，但是留了个悬念；生物信息中第一个要掌握的基础技能，没错，就是这一讲的文件操作了。
本文就简单介绍一下，生物信息中那些基础的文件读写，后续补充文件权限管理以及文件的移动、复制等操作（先画个饼，有空在填坑）。
文件读写
在日常的生物信息学分析中，我们一般都是直接面对一些文件进行操作的，比如测序数据 fastq文件，在比如比对结果bam或者sam文件，还有 gtf或者gff3格式的注释文件甚至还有一些数据库文件等等。
面对这些文件时候交互式编程会很力不从心，大多数情况下我们都会选择脚本式编程。（关于交互式编程与脚本式编程的区别见另一篇文章，当然现在还没有写！(lll￢ω￢)）。
脚本式编程中首先遇到的就是文件的读写，其次是参数的传递。
一般的文件都可以通过open的方式打开，readline的方式进行读取。
# 读文件
f = open('/path/to/file','r') # 只读模式打开
line = f.readline() # 读取一行
linse = f.readlines() # 读取所有行
f.close() # 关闭打开的文件
# 写文件
f = open('/path/to/file','w') # 写模式打开
f.write(string) # 将 string 写入文件中
f.close() # 关闭打开的文件
这里推荐使用with语句,不光可以减少代码量，还可以处理一些异常情况：
with open('/path/to/file','r') as f:
line = f.readline() # 读取一行
linse = f.readlines() # 读取所有行
生物信息中一般主流的文件格式有fasta、fastq、sra、Genebank、gtf、gff3、sam、bam等。
fasta文件读写
fasta文件一般以 > 开头，第一行一般为序列名称，或者其他信息之类的，第二行开始为其碱基序列。如下：
# test.fa
>MSTR.8.1 gene=MSTR.8
CCACAATTAATCAGAGGTAGCGACTACCCTATTTCTTCATTTCTTATGATTCTTACTATTTAAGCACTCA
TTTTCCATACCCCCACCTTTTACTTTCCTCTATTTTCCTCTCATCGTTAACTGATTCTGTTCAAGTTATA
GGGTTGGTTCTCAATCCTTGATATGTTGCTGGAATCCCTTTGGATTCGCTAAAAGGAGGATTTCAAGGGT
TTTCCAGAGTTTCCGACATTTGATTAAGGCGATCGGTGGATTTGAACTTGATTTCGATTGTTTCACAGTA
CCACCATGTATGGGTAATCGAGGATTTTAGGGAGGAGAGATTTGAGCCTCAAGGAGGAAGATAACCAGTT
GGTATTTCATGGTAGTAGAATCTTCAGCAATTTCTAGGTATGTTTAGAGGGT
>MSTR.26.1 gene=MSTR.26
CACCAGATCAAGATGACAGTAACTACAATTTAAGCAGATCGAACAGGAAACATGTAGCCATCTTTCATGT
ATGAAAGGTCCAAATACTAGACTCTACTTGTGCATCCCTCAGACTACATTTATTTATCTTGCCGAAATGT
TAAAAATGAAGTTTCAGTTTCTCATGTAAAGTTTAGTTTCAGTTTCATTGAGATGAATCCATGGATGTAG
CTAGTCATCCATTCTCGGGAAACAACAACCACC
对于fasta文件的读取方式，一般有以下两种方法：
- 直接读取
将fasta文件的每一条序列以字典的形式储存起来，字典的键为fasta序列的名称，值为其碱基序列。
#!/usr/bin/evn python
# -*- coding: UTF-8 -*-
with open("test.fa",'r') as fa:
sequences = {}
for line in fa:
if line.startswith(">"):
name = line.rstrip("\n")
sequences[name] = ""
else:
sequences[name] = sequences[name] + line.rstrip("\n")
print (sequences)
输出如下：
{'>MSTR.26.1 gene=MSTR.26': 'CACCAGATCAAGATGACAGTAACTACAATTTAAGCAGATCGAACAGGAAACATGTAGCCATCTTTCATGTATGAAAGGTCCAAATACTAGACTCTACTTGTGCATCCCTCAGACTACATTTATTTATCTTGCCGAAATGTTAAAAATGAAGTTTCAGTTTCTCATGTAAAGTTTAGTTTCAGTTTCATTGAGATGAATCCATGGATGTAGCTAGTCATCCATTCTCGGGAAACAACAACCACC', '>MSTR.8.1 gene=MSTR.8': 'CCACAATTAATCAGAGGTAGCGACTACCCTATTTCTTCATTTCTTATGATTCTTACTATTTAAGCACTCATTTTCCATACCCCCACCTTTTACTTTCCTCTATTTTCCTCTCATCGTTAACTGATTCTGTTCAAGTTATAGGGTTGGTTCTCAATCCTTGATATGTTGCTGGAATCCCTTTGGATTCGCTAAAAGGAGGATTTCAAGGGTTTTCCAGAGTTTCCGACATTTGATTAAGGCGATCGGTGGATTTGAACTTGATTTCGATTGTTTCACAGTACCACCATGTATGGGTAATCGAGGATTTTAGGGAGGAGAGATTTGAGCCTCAAGGAGGAAGATAACCAGTTGGTATTTCATGGTAGTAGAATCTTCAGCAATTTCTAGGTATGTTTAGAGGGT'}
- 使用第三方模块
使用 Biopython SeqIO模块，不多使用前好像需要先安装一下，怎么安装请自行google。
#!/usr/bin/evn python
# -*- coding: UTF-8 -*-
from Bio import SeqIO
with open ("test.fa",'r') as fa:
for record in SeqIO.parse(fa,'fasta'):
print(record)
输出如下：
ID: MSTR.8.1
Name: MSTR.8.1
Description: MSTR.8.1 gene=MSTR.8
Number of features: 0
Seq('CCACAATTAATCAGAGGTAGCGACTACCCTATTTCTTCATTTCTTATGATTCTT...GGT', SingleLetterAlphabet())
ID: MSTR.26.1
Name: MSTR.26.1
Description: MSTR.26.1 gene=MSTR.26
Number of features: 0
Seq('CACCAGATCAAGATGACAGTAACTACAATTTAAGCAGATCGAACAGGAAACATG...ACC', SingleLetterAlphabet())
正确读取后就可以进行其他需求的操作了。
fastq文件读写
fastq文件格式：，每条序列有4行，第一行一般以 @ 符号开头，主要是测序的一些信息可以简单理解为序列名称，第二行是其碱基序列，第三行也是名称，一般为了节省空间就用 + 号表示，第四行是序列每个碱基的测序质量值，与第二行的碱基一一对应。
@FCD1VGJACXX:2:1101:5762:2000#AAGTCGCG/1
NGGGAGACCGGGGTTCGATTCCCCGACGGGTCGTATGCCGTCTTCTGCT
+
BS\ceeeegggggffgighiiiihhUegiigggeeeeddcccccccccc
@FCD1VGJACXX:2:1101:7844:2000#AAGTCGCG/1
NTGGCTTCCTAAGCCAGGGATTGTGGGTTCGTATGCCGTCTTCTGCTTG
+
BPYa^accgeeggihfhhhfhiQQQQPPPPPP^XYefgeeWaeghhgih
对于fastq文件我们大部分情况下主要考虑其碱基质量，常规情况下，会对其序列进行碱基统计、质量统计等，对于fastq文件，我们同样也有两种方式可以进行读取：
1. 常规方式
#!/usr/bin/evn python
# -*- coding: UTF-8 -*-
# style 1
with open("test.fq",'r') as fq:
lines = fq.readlines()
head = [item[:-1] for item in lines[::4]]
read = [item[:-1] for item in lines[1::4]]
qual = [item[:-1] for item in lines[3::4]]
sequence = dict(zip(read, qual))
print(sequence)
# style 2
with open('test.fq','r') as fq:
sequence = {}
name = None
for i,line in enumerate(fq):
if i % 4 == 0:
name = line.rstrip('\n')
if i % 4 == 1:
read = line.rstrip('\n')
if i % 4 == 3:
sequence[(name,read)] = line.rstrip('\n')
for key in sequence:
print (key[0],key[1],sequence[key])
输出如下：
# style 1
{'NTGGCTTCCTAAGCCAGGGATTGTGGGTTCGTATGCCGTCTTCTGCTTG': 'BPYa^accgeeggihfhhhfhiQQQQPPPPPP^XYefgeeWaeghhgih', 'NGGGAGACCGGGGTTCGATTCCCCGACGGGTCGTATGCCGTCTTCTGCT': 'BS\\ceeeegggggffgighiiiihhUegiigggeeeeddcccccccccc'}
# style 2
('@FCD1VGJACXX:2:1101:5762:2000#AAGTCGCG/1', 'NGGGAGACCGGGGTTCGATTCCCCGACGGGTCGTATGCCGTCTTCTGCT', 'BS\\ceeeegggggffgighiiiihhUegiigggeeeeddcccccccccc')
('@FCD1VGJACXX:2:1101:7844:2000#AAGTCGCG/1', 'NTGGCTTCCTAAGCCAGGGATTGTGGGTTCGTATGCCGTCTTCTGCTTG', 'BPYa^accgeeggihfhhhfhiQQQQPPPPPP^XYefgeeWaeghhgih')
2. 第三方模块读取
与读取fasta文件的模块一样，使用 Biopython
#!/usr/bin/evn python
# -*- coding: UTF-8 -*-
from Bio import SeqIO
with open ("test.fq",'r') as fq:
for record in SeqIO.parse(fq,'fastq'):
print(record)
输出如下：
ID: FCD1VGJACXX:2:1101:5762:2000#AAGTCGCG/1
Name: FCD1VGJACXX:2:1101:5762:2000#AAGTCGCG/1
Description: FCD1VGJACXX:2:1101:5762:2000#AAGTCGCG/1
Number of features: 0
Per letter annotation for: phred_quality
Seq('NGGGAGACCGGGGTTCGATTCCCCGACGGGTCGTATGCCGTCTTCTGCT', SingleLetterAlphabet())
ID: FCD1VGJACXX:2:1101:7844:2000#AAGTCGCG/1
Name: FCD1VGJACXX:2:1101:7844:2000#AAGTCGCG/1
Description: FCD1VGJACXX:2:1101:7844:2000#AAGTCGCG/1
Number of features: 0
Per letter annotation for: phred_quality
Seq('NTGGCTTCCTAAGCCAGGGATTGTGGGTTCGTATGCCGTCTTCTGCTTG', SingleLetterAlphabet())
至于sam/bam文件的读取,可以使用第三方模块：pysam ，具体练习在后面遇到了在说。
gtf及gff3注释信息文件的读取与一般文本文件的读取方式一致。
