1 python生信入门

1. 解决的实际问题，计算人类参考基因组hg38每条染色体碱基数据ATCG的比例

2. 计算人类参考基因组hg38的外显子区域碱基数量

3. 目的是了解基础python语法

4. 参考书目：python编程从入门到实践，参考视频资料见最底部

5. 测试数据下载地址链接：https://share.weiyun.com/5Xmmzdz 密码：6q2y5j

打开python3

vim 1.txt
:wq  ##Esc+ ：wq保存退出
python3 # 打开python3

f = open('./1.txt','r')
for i in f:
  print (i)

with open('./1.txt') as file_object:
  for line in file_object:
   print (line)

图1

f = open('./2.txt','w')

with open('./1.txt','r') as imput_f:
  for line in imput_f:
   line = line.strip()
   result= int(line)+int(line)
   f.write(str(result)+'\n')

f.close()

图2

output_file = open('./2.txt','a')  ## 附加内容

with open('./1.txt','r') as imput_f:
  for line in imput_f:
   line = line.strip()
   result= int(line)+int(line)
   output_file.write(str(result)+'\n')


output_file.close()

图3

完成探索人类基因组GC含量与gene密度的关系；

提示：

1.将染色体成分1Mb大小的bin统计GC含量；

2.将染色体分为1Mb大小的bin统计gene的密度；

3.使用R将这两者用散点图进行绘图；

4.使用R的线性回归函数lm()进行回归；

5.R回归部分，参考资料

blog.fens.me/r-linear-regression/

https://github.com/menghaowei

def my_abs(imput):
     if imput >= 0:
         return imput
     else:
         return -imput
        
my_abs(-10)

作业

人类的线粒体长度大约为16.7Kb，请读取人类线粒体参考基因组序列信息（fasta格式），统计各种碱基的数量以及人类线粒体参考基因组的总长度。

编程作业

# _*_ coding:UTF-8 _*_
### 保证中文不报错pycharm

选取参考基因组hg38作为测试数据1000-1300行 n 和 p是同时出现的
cat ~/reference/UCSC/hg38/hg38.fa| sed -n '1000,1300p' >> ~/project/python/t.txt
cd - # 回到上一个目录
python3

图4

chrM_genome = ''
with open('./t.txt','r') as imput_f:
  for line in imput_f:
   line = line.strip().upper()  ## 去除前后空格，大写
   chrM_genome = chrM_genome + line
    
print (len(chrM_genome))
print (chrM_genome.count('A'))
print ('the char A count is : {0}'.format(chrM_genome.count('A')))  # format 将其格式化

Fastq 格式 四行

图5

>>> with open(r'./fa_q.txt','r') as imput:
...  for index,line in enumerate(imput):
...   if index % 4==0:
...    print ('>'+line.strip()[1:])
...   elif index % 4 ==1:
...    print (line.strip())
...   elif index % 4==2:
...    continue
...   elif index %4==3:
...    continue
... 
>SRR1039512.325 HWI-ST177:314:C12K6ACXX:1:1101:5657:2438 length=63
TCAGAAGAAAAGCATATGTCTGCCATGGGACTATTGCAGTGCGTCTCCATCAGTGTTAACACA
>SRR1039512.326 HWI-ST177:314:C12K6ACXX:1:1101:5550:2446 length=63
TAAGCAGTGCTTGAATTATTTGGTTTCGGTTGTTTTCTATTAGACTATGGTGAGCTCAGGTGA
>SRR1039512.327 HWI-ST177:314:C12K6ACXX:1:1101:5908:2367 length=63
TGGGGAATAGCTTCTCAAGGGCAAGTGCTGTTTCCATACTAGTTCTTTTCCTTGCTCTCTTAT
>SRR1039512.328 HWI-ST177:314:C12K6ACXX:1:1101:5852:2396 length=63
GGAGGAAGCCCACGGTGGGATCTGGCACAAAGATGTAGAGCCGTTTCCGCTCATCGGTCTCCA
>SRR1039512.329 HWI-ST177:314:C12K6ACXX:1:1101:5768:2409 length=63
ACCTGGCAGTGCTGCAACAGTATCTGCCTCTACAAGTAACATCATACCCCCAAGACACCAGAA
>SRR1039512.330 HWI-ST177:314:C12K6ACXX:1:1101:5840:2413 length=63
GCAGAGGTAATTGCCTGAATCCTTCTGTTGTAGACTACGTAGCAGAAGGCCTTGATCTGTCCT
>SRR1039512.331 HWI-ST177:314:C12K6ACXX:1:1101:5826:2463 length=63
GAAAGGCAAGGAGGAAGCTTATCTATGAAAAAGCAAAGCACTATCACAAGGAATATAGGC
>SRR1039512.332 HWI-ST177:314:C12K6ACXX:1:1101:5888:2477 length=63
GTCTGTCTCTTGCTTCAACAGTGTTTGGACGGAACAGATCCGGGGACTCTCTTCCAGCCTCCG
>SRR1039512.333 HWI-ST177:314:C12K6ACXX:1:1101:6208:2288 length=63
GTCATAACCACCACCGCTTACACCTGGAGGTCCAGGAGGGCCAGGGGGACCAGGGGGGCCAGC

##每隔40一行
with open(r'C:\Users\kexin\Desktop\python\fa_q.txt','r') as imput:
  for index,line in enumerate(imput):
   if index % 4==0:
    print ('>'+line.strip()[1:])
   elif index % 4 ==1:
    for i in range(0,len(line.strip()),40):
     print (line.strip()[i:i+40])
   elif index % 4==2:
    continue
   elif index %4==3:
    continue

output_file = open (r'C:\Users\kexin\Desktop\python\result2.fastq','w')
with open(r'C:\Users\kexin\Desktop\python\fa_q.txt','r') as imput:
  for index,line in enumerate(imput):
   if index % 4==0:
    output_file.write(('>'+line.strip()[1:])+'\n')
   elif index % 4 ==1:
    for i in range(0,len(line.strip()),40):
     output_file.write(line.strip()[i:i+40]+'\n')
   elif index % 4==2:
    print (index)
   elif index %4==3:
    continue
output_file.close()

读取fasta文件

with open (r'C:\Users\kexin\Desktop\python\result2.fastq','r') as input_file:
    seq=''
    header= input_file.readline().strip()[1:]  ## 读取第一行
    for line in input_file:
        if line[0]!='>':
            seq =seq.strip()+line  ### 这里是要确保刚刚被40个一行分隔的seq数据重新叠加
        elif line[0]=='>':
            print (header)
            print (seq)
            header = line.strip()[1:]
            seq = ''       ### 一个循环结束将seq清空
    print (header)
    print (seq)

一个能读取fasta格式的函数（收藏版）

def read_fasta(file_path=''):
    line = ''
    try:
        fasta_handle=open(file_path,'r')
    except:
        raise IOError('you input FAST file is not right')
    while True:
        line = fasta_handle.readline()
        if line == '':
            return
        if line[0] =='>':
            break
 ## while the file is not empty, load FASTA file   
    while True:
        if line[0] !='>':
            raise ValueError("Records in Fasta files should start with '>'")
        title = line[1:].rstrip()
        lines =[]
        line = fasta_handle.readline()
        while True:
            if not line:
                break
            if line[0] == '>':
                break
            lines.append(line.rstrip())
            line = fasta_handle.readline()
        yield title, ''.join(lines).replace('','').replace('\r','')
        if not line:
            return
    fasta_handle.close()
    assert False, 'your input fasta file have format problem'
    
    
read_fasta(file_path=r'C:\Users\kexin\Desktop\python\result2.fastq')  ## 出现一个报错，在windows里面，要加上r，不然\代表转义:参考https://blog.csdn.net/caibaoH/article/details/78335094

for header,seq in read_fasta(file_path=r'C:\Users\kexin\Desktop\python\result2.fastq'):
    print (header)
    print (seq)

计算人类参考基因组的序列长度

分析步骤：

利用读取fasta的代码读取基因组序列信息

统计每条序列的长度，并输出

统计人类参考基因组的总长度，并输出

统计每条染色体N的长度，并输出

提取基因组特定区域的序列，并输出（支持+/-链）

hg38_genome={}
for chr_name, chr_seq in read_fasta(file_path=r'/four/yqchen/reference/UCSC/hg38/hg38.fa'):
    hg38_genome[chr_name]=chr_seq

这里要用绝对路径不然报错

图6

hg38_genome.keys()
for chr_name in hg38_genome.keys():
    print (chr_name,len(hg38_genome.get(chr_name)))
## 这里有点不对，但是无所谓了

图7

hg38_total_len=0
for index in range (1,26):
    if index<= 22:
        chr_name = 'chr{0}'.format(index)
    elif index == 23:
        chr_name = 'chrX'
    elif index == 24:
        chr_name = 'chrY'
    elif index == 25:
        chr_name = 'chrM'
    print (chr_name,len(hg38_genome.get(chr_name)))
    

这里结果是正确的，这是为什么？这里没搞明白的。

图8

hg38_total_len = 0
for index in range (1,26):
    if index<= 22:
        chr_name = 'chr{0}'.format(index)
    elif index == 23:
        chr_name = 'chrX'
    elif index == 24:
        chr_name = 'chrY'
    elif index == 25:
        chr_name = 'chrM'
    hg38_total_len= hg38_total_len + int(len(hg38_genome.get(chr_name))) 

print (hg38_total_len) ### 和上面那一句分开执行

图9
图10

len(hg38_genome['chr21'])
hg38_genome['chr21'].count('N')
hg38_genome['chr22'].count('N')

## 正负链

chr_name='chr1'
chr_start=10000
chr_end=10300
hg38_genome[chr_name][chr_start-1:chr_end-1].upper()

#### 反向互补
import string
translate_table=string.maketrans('ATGC','TACG')
a=hg38_genome[chr_name][chr_start-1:chr_end-1].upper()
a.translate(translate_table)[::-1]

##方向互补定义函数
import string
def com_rev(input):
    translate_table=string.maketrans('ATGC','TACG')
    return input.upper().translate(translate_table)[::-1]

编码基因长度，从参考基因组注释信息找

with open(r'C:/Users/kexin/Desktop/python/hg38_RefSeq_genetable') as input_f:
    header = input_f.readline()
    for line in input_f:
        line=line.strip().split('\t')
        print (line)
        break

with open(r'C:/Users/kexin/Desktop/python/hg38_RefSeq_genetable') as input_f:
    header = input_f.readline()   #### 读一行 不要用readlines
    for line in input_f:
        line=line.strip().split('\t')
        exon_count = int(line [8])
        exon_start =line[9].strip(',') ### 去掉最后的逗号
        print (exon_start)
        break

with open(r'C:/Users/kexin/Desktop/python/hg38_RefSeq_genetable') as input_f:
    header = input_f.readline()   #### 读一行 不要用readlines
    for line in input_f:
        line=line.strip().split('\t')
        # print (line)
        exon_count = int(line [8])
        exon_start =line[9].strip(',')
        exon_end =line[10].strip(',')### 去掉最后的逗号
        print (exon_start)
        print (exon_end)
        break

with open(r'C:/Users/kexin/Desktop/python/hg38_RefSeq_genetable') as input_f:
    header = input_f.readline()   #### 读一行 不要用readlines
    for line in input_f:
        line=line.strip().split('\t')
        # print (line)
        exon_count = int(line [8])
        exon_start = line[9].strip(',').split(',') ### map函数int用于后面的列表
        exon_end = line[10].strip(',').split(',')
        print (exon_start)
        print (exon_end)
        break

with open(r'C:/Users/kexin/Desktop/python/hg38_RefSeq_genetable') as input_f:
    header = input_f.readline()   #### 读一行 不要用readlines
    for line in input_f:
        line=line.strip().split('\t')
        # print (line)
        exon_count = int(line [8])
        exon_start = map(int,line[9].strip(',').split(',')) ### map函数int用于后面的列表
        exon_end = map(int,line[10].strip(',').split(','))
        print (exon_start)
        print (exon_end)
        break

with open(r'C:/Users/kexin/Desktop/python/hg38_RefSeq_genetable') as input_f:
    header = input_f.readline()   #### 读一行 不要用readlines
    for line in input_f:
        line=line.strip().split('\t')
        # print (line)
        exon_total_len=0
        exon_count = int(line [8])
        exon_start = map(int,line[9].strip(',').split(',')) ### map函数int用于后面的列表
        exon_end = map(int,line[10].strip(',').split(','))
        
        for index in range(exon_count):
            exon_total_len= exon_total_len+exon_end[index]-exon_start[index]
        print (exon_total_len)
        print (exon_start)
        print (exon_end)
        break

这里python3 会报错

1555423217495.png

解决方案，在map基础上加上list

https://blog.csdn.net/mingyuli/article/details/81238858

1555423298540.png

with open(r'C:/Users/kexin/Desktop/python/hg38_RefSeq_genetable') as input_f:
    header = input_f.readline()   #### 读一行 不要用readlines
    for line in input_f:
        line=line.strip().split('\t')
        # print (line)
        exon_total_len=0
        exon_count = int(line [8])
        exon_start = list(map(int, line[9].strip(',').split(','))) ### map函数int用于后面的列表
        exon_end = list(map(int, line[10].strip(',').split(',')))
        
        for index in range(exon_count):
            exon_total_len= exon_total_len+exon_end[index]-exon_start[index]
        print (exon_total_len)
        print (exon_start)
        print (exon_end)
        break

完整计算版：

gene_exon_len_dict={}
with open(r'C:/Users/kexin/Desktop/python/hg38_RefSeq_genetable') as input_f:
    header = input_f.readline()   #### 读一行 不要用readlines
    for line in input_f:
        line=line.strip().split('\t')
        # print (line)
        exon_total_len=0
        exon_count = int(line [8])
        exon_start = list(map(int, line[9].strip(',').split(','))) ### map函数int用于后面的列表
        exon_end = list(map(int, line[10].strip(',').split(',')))
         
        for index in range(exon_count):
            exon_total_len= exon_total_len+exon_end[index]-exon_start[index]
        gene_id = line[12]
        
        if gene_exon_len_dict.get(gene_id) is None: #如果gene_exon_len_dict字典里没有这个基因，加入这个基因的长度
            gene_exon_len_dict[gene_id] = exon_total_len
        else:
            if gene_exon_len_dict.get(gene_id) < exon_total_len:
                gene_exon_len_dict[gene_id] = exon_total_len   

len(gene_exon_len_dict)

exon_total =0 
for gene_id in gene_exon_len_dict.keys():
    exon_total = exon_total+gene_exon_len_dict[gene_id]
    # exon_total = exon_total+gene_exon_len_dict.values(gene_id) 这是错误写法
    
print (exon_total)


exon_rate = 111336157 / 3088286401 *100
exon_rate

推荐视频

https://www.bilibili.com/video/av16064125?from=search&seid=5895915127692784448

https://www.bilibili.com/video/av15887563?from=search&seid=5895915127692784448