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

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

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

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

打开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/

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

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

推荐视频