在学习了生信大神孟浩巍的知乎Live “学习Python, 做生信”之后,对第二部分的文件信息处理部分整理了如下的笔记。
一、fasta与fastq格式的转换
1、首先需要了解FASTA和FASTQ格式的详解
1)具体的详解看知乎专栏的这篇文章,写的很详细。
https://zhuanlan.zhihu.com/p/20714540
2)关于FASTA
- 主要分为两部分:第一行是“>”开始的储运存的序列描述信息;接下来是序列部分。序列部分可以有空格,按照一般70~80bp。用python进行操作的过程中,需要去掉空格和换行符,把序列读成一行再处理。(.strip()可以除去空格)
- 即使是mRNA的序列,为了保证数据的统一性,序列中的U依然是用T来表示。核苷酸序列信息对应列表如下所示。
A adenosine C cytidine G guanine
T thymidine N A/G/C/T (any) U uridine
K G/T (keto) S G/C (strong) Y T/C (pyrimidine)
M A/C (amino) W A/T (weak) R G/A (purine)
B G/T/C D G/A/T H A/C/T
V G/C/A - gap of indeterminate length
3)关于FASTQ
- 刚刚测到的测序数据一般用FASTQ格式进行储存,因为其中包含了测序质量信息。
- 四行依次为:信息头,序列信息,注释信息,质量信息。质量值的计算方法见上面链接里的文章。
- 序列符号为 ATCGN,N表示无法确定。
2、读取fastq文件并进行相应的格式转化
1)主要流程如下:
- 读取fastq文件,第1行的header作为fasta的header
- 读取fastq文件,第2行的seq作为fasta格式的seq
- 第3/4行直接忽略
- 格式化输出fasta文件
2)以下是没有设置输出并储存为.fa文件的代码
1 #-*- coding: UTF-8 -*-
2 with open( 'E:\genome_file\\test.fastq','r') as input_fastq:
3 for index, line in enumerate (input_fastq):
4 if index % 4 == 0:
5 print ">" + line.strip()[1:]
6 elif index % 4 == 1:
7 for i in range (0,len(line.strip()),40):
8 print line.strip()[i:i+40]
9 elif index % 4 == 2:
10 pass
11 elif index % 4 == 3:
12 pass
3)以下是设置了输出fasta文件的代码:
1 out_file = open(E:\genome_file\\test.fa)
2
3 with open( 'E:\genome_file\\test.fastq','r') as input_fastq:
4
5 for index, line in enumerate (input_fastq):
6
7 if index % 4 == 0:
8 out_file.write(">" + line.strip()[1:]+"\n" )
9
10 elif index % 4 == 1:
11 for i in range (0,len(line.strip()),40):
12 out_file.write(line.strip()[i:i + 40] +"\n")
13
14 elif index % 4 == 2:
15 pass
16
17 elif index % 4 == 3:
18 pass
4)有几点注意的地方这里提一下:
- 2)和3)中代码的主要区别就是在3)在一开始创建了一个output_file,这就是用来输出用的。所以可以看到2)中的输出直接用的是print,而3)中则是通过.write写入到 output_file当中。
- 读取第一行的时候要特别注意:去掉原先FASTQ文件中的“@”字符,并加上一个“>”字符。
- enumarate是个遍历函数,能够输出相应的索引和对应的值。
- print能够自动换行,但是write不行。所以在写入的时候需要切片+ \n,切多少呢一般?前面在关于fasta格式文件介绍中有提到。
二、读取fasta格式文件
1、分析过程
- 当有>的时候,就为标题行;
- 当不是>的时候,就是序列信息;
- 当是序列信息的时候,需要进行序列的拼接;
- 最终返回序列的header与seq
2、简化版的主要代码如下:
input_file = open( 'E:\genome_file\\test.fa','r')
seq = ""
header = input_file.readline().strip()[1:]
for line in input_file:
line = line.strip()
if line[0] != ">":
seq = seq + line
else:
print header
print seq
header = line[1:]
seq = ""
#打印最后一行
print header
print seq
input_file.close()
3、值得注意的几个地方:
- 整体思路就是通过for循环进行序列信息的累加;
- 一开始先读一行header(虽然暂时也还没搞太明白为什么,反正就是如果不读的话出来结果是乱七八糟的)
- 按照上面那个情况,循环中是无法打印最后一行的,所以要在最外面打印一下。
- readline和readlines的区别:readline只读一行,readlines则是一下子读整个文件。永远不要用readlines!!
4、下面的代码是汇总了上面的功能,直接通过函数形式fastq-->fasta格式文件的转换
1 def read_fasta(file_path=""):
2 """
3 Loading FASTA file and return a iterative object
4 """
5
6 line = ""
7
8 try:
9 fasta_handle = open(file_path,"r")
10 except:
11 raise IOError("Your input FASTA file is not right!")
12
13 # make sure the file is not empty
14 while True:
15 line = fasta_handle.readline()
16 if line == "":
17 return
18 if line[0] == ">":
19 break
20
21 # when the file is not empty, we try to load FASTA file
22 while True:
23 if line[0] != ">":
24 raise ValueError("Records in Fasta files should start with '>' character")
25 title = line[1:].rstrip()
26 lines = []
27 line = fasta_handle.readline()#这里是读第二行了
28 while True:#循环只用来加载序列行
29 if not line:
30 break
31 if line[0] == ">":
32 break
33 lines.append(line.rstrip())
34 line = fasta_handle.readline()
35
36 yield title,"".join(lines).replace(" ","").replace("\r","")
37
38 if not line:
39 return
40
41 fasta_handle.close()
42 assert False, "Your input FASTA file have format problem."
43
44 for header,seq in read_fasta(file_path=r"D:\data_file\test.fa"):
45 print header
46 print seq
47 #下面是后边的练习中以hg19为例进行操作
48 hg19_genome = {}
49 for chr_name , chr_seq in read_fasta(file_path=r"D:/data_file/hg19_only_chromosome.fa"):
50 hg19_genome[chr_name] = chr_seq
其实不太明白的是这里边最终是返回两个head和seq两个值了吗?为什么后边直接来两个参数就能开始for循环?
三、计算人类基因组长度信息
1、下载人类参考基因组信息(UCSC网站)
- http://hgdownload.soe.ucsc.edu/downloads.html#human
- 下载对象为hg19全基因组信息
2、利用读取fasta的代码读取基因组序列
1)代码如下:
1 #前面已经定义过的read_fasta函数这里不再重复写。
2 hg19_genome = {}
3
4 for chr_name , chr_seq in read_fasta(file_path=r"D:/data_file/hg19_only_chromosome.fa"):
5 hg19_genome[chr_name] = chr_seq
2)注意几点
- 程序中把下载的.fa文件中的信息输入到hg19_genome的列表当中
- 读取一个全基因组数据需要耗费一定量的时间和内存,如果再pycharm或者通过powershell进行运行,调试的时候每调试一次就要运行一次,很费事。这个时候jupyter就方便多了。使用方法:1、终端输入conda,回车;2、输入jupyter notebook,回车;3、在弹出的中选择计算机上相应的Python程序,hello world。
3、统计每条序列的长度,并输出 ;
1)下面这个是一个乱序的输出
1 for chr_name in sorted(hg19_genome.keys()):
2 print chr_name, len(hg19_genome.get(chr_name))
2)下面的代码能够做到按顺序输出
1 hg19_total_len = 0
2 for index in range (1,26):
3 if index <=22:
4 chr_name = "chr%s" %index
5 elif index == 23:
6 chr_name = "chrX"
7 elif index == 24:
8 chr_name = "chrY"
9 elif index == 25:
10 chr_name = "chrM"
11 print chr_name, len(hg19_genome.get(chr_name))
12 hg19_total_len = hg19_total_len + len(hg19_genome.get(chr_name))
13
14 print "Total chromosome length is %s" %hg19_total_len
注意一点:一般从字典里边根据Key来获取内容,用.get(),这样子在对象不存在的时候就不会报错。
3、统计人类参考基因组的总长度,并输出
在上面已经输出,往回看。
4、统计每条染色体N的长度,并输出(GC含量同理)
以N为例,其他同理。使用.count指令。
1 hg19_genome['chr22'].count("N")
5、提取基因组特定区域的序列,并输出(支持+/-链)
1)思路分析
- 通过切片可以截取特定区域的序列
- 需要进行转移来过得互补序列
- [::-1]能够实现序列反向的功能
2)上代码
1 chr_name = "chr1"
2 chr_start = 10000
3 chr_end = 10100
4 #特定区域截取
5 hg19_genome[chr_name][chr_start-1:chr_end-1].upper()
6 #互补链
7 import string
8 translate_table = string.maketrans("ATCG","TAGC")
9 a = hg19_genome[chr_name][chr_start-1:chr_end-1].upper()
10 a.translate(translate_table)
11
12 #反向链
13 a.translate(translate_table)[::-1]#反向
- 这里导入string模块,设定一个转译表
- .upper是用来把小写字母全部转换成大写
- 注意要“-1”
3)定义成函数方便使用
1 def con_rev(input_seq):
2 translate_table = string.maketrans("ATCG","TAGC")
3 output_seq = input_seq.upper().translate(translate_table)[::-1]
4 return output_seq
5
6 con_rev("aaTTCCGG")
四、计算人类基因组外显子的长度(CDS区域同理)
1、下载人类参考转录组(UCSC table browser)
- 一般在track一栏均选择RefSeq Genes
2、整体思路
- 从转录组抓取关键信息:exon_count, gene_id, exon_start,exon_end
- 最后输出一个字典:{‘基因名,外显子长度’}
- 同一基因存在不同转录本,本例中按照最长计算,所以需要进行比较
3、上代码
1 gene_exon_len_dict = {}
2 with open(r"E:\genome_file\hg19_refseq_gene_table","r") as input_file:
3 header = input_file.readline()
4 for line in input_file:
5 line = line.strip().split("\t")
6
7 exon_count = int(line[8])
8 exon_start = line[9].split(",")#还有更好的方法进行index
9 exon_end = line[10].split(",")
10 exon_total_len = 0
11 gene_id = line[12]
12 for index in range (exon_count):
13 exon_total_len = exon_total_len + int(exon_end[index]) - int(exon_start[index])
14
15 if gene_exon_len_dict.get(gene_id) is None:#注意:如果直接用dict[gene_id]是会出现key error的
16 gene_exon_len_dict[gene_id] = exon_total_len
17
18 elif gene_exon_len_dict[gene_id] < exon_total_len:
19 gene_exon_len_dict[gene_id] = exon_total_len
20
21 else:
22 pass
23
24 #注意啊调试的时候要加break,这是很重要的!!
25
26 print gene_exon_len_dict
27
- .split()能够根据括号中的内容把相应的字符串拆分成为列表
- exon_start和exon_end都要以整型的格式才能进行相减。除了上面代码中的处理方法,还可以通过如下方法来实现:
1 exon_start = map(int,line[9].strip(",").split(","))
2 exon_end = map(int,line[10].strip(",").split(","))
- 特别重要!!数据量大的运算,调试的时候,在后边加个 break啊!!不然卡的你怀疑人生~~~