在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