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

在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