人类基因组本地化及简单分析

在NCBI上下载 GRCh38

wget ftp://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/annotation/GRCh38_latest/refseq_identifiers/GRCh38_latest_genomic.fna.gz

解压文件（.fasta, .fa, .fna, .fsa, .mpfa）

gzip
-d GRCh38_latest_genomic.fna.gz
#人的h38基因组是3G的大小，一个英文字符是一个字节，所以30亿bp的碱基就是3G左右

head GRCh38_latest_genomic.fna

查看该文件可以看到，里面有很多的N，这是基因组里面未知的序列，用N占位，但是觉得部分都是A.T.C.G这样的字符，大小写都有，分别代表不同的意思

统计了一下里面这个文件的行数

time wc -l GRCh38_latest_genomic.fna

用awk统计行数（效率相比wc –l 慢）

time awk 'END { print NR }' GRCh38_latest_genomic.fna

看一下标题行

grep '>' GRCh38_latest_genomic.fna | sed -n 'p'
grep '>' GRCh38_latest_genomic.fna | sed -n 'p' >> list.txt

统计每个标题下基因片段的长度，提取标题和长度写入一个新文件

time python GECh38_title_length.py

fasta_file=open('/home/sunchengquan/GRCh38_latest_genomic.fna','r')
out_file = open('GRCh38_title_length.txt','w')
seq = ''
i = 0
for line in fasta_file:
if line[0] == '>' and seq == '':
header = line.strip()
elif line[0] != '>':
seq =seq + line.strip()
elif line[0] == '>' and seq != '':
num = len(seq)
out_file.write(header +'n'+ str(num)+ 'n')
i += 1
print('writing:',i)
seq = ''
header = line.strip()
out_file.close()

看一下GRCh38_title_length.txt里面的内容

提取标题行，添加到列表，并打印

time python GECh38_title.py

input_file=open("/home/sunchengquan/GRCh38_latest_genomic.fna","r")
title_list = []
for line in input_file:
if line[0] == '>':
field = line
title_list.append(field)
print(field)
类似于
grep '>' GRCh38_latest_genomic.fna | sed -n 'p' > list.txt

最后

以上就是热情橘子为你收集整理的人类基因组本地化及简单分析的全部内容，希望文章能够帮你解决人类基因组本地化及简单分析所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。