Python 入门之文件操作(一)

上一讲中我们初步认识了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文件的读取方式,一般有以下两种方法:

  1. 直接读取
    将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'}
  1. 使用第三方模块
    使用 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注释信息文件的读取与一般文本文件的读取方式一致。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 206,214评论 6 481
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 88,307评论 2 382
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 152,543评论 0 341
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 55,221评论 1 279
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 64,224评论 5 371
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,007评论 1 284
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,313评论 3 399
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,956评论 0 259
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 43,441评论 1 300
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,925评论 2 323
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,018评论 1 333
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,685评论 4 322
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,234评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,240评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,464评论 1 261
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,467评论 2 352
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,762评论 2 345

推荐阅读更多精彩内容