Bioinfo:学习Python，做生信PartII 学习笔记

Posted 2020-10-19 dingtou00

在学习了生信大神孟浩巍的知乎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啊！！不然卡的你怀疑人生~~~