在进行单细胞测序的时候,我们通常最后得到的是 cell*geneid 的 matrix 文件,有时候我们需要 geneid 以及 geneid 对应的 fasta 序列数据,这就需要让我们做一些处理。

我们所需的文件是 gtf 文件和 基因组文件,需要使用的工具是 gffread 文件,这些大家都可以从官网进行下载。

以我们的数据为例,我们使用的是 GRCh38.p13.genome.fa 的数据,所以便是可以进行如下的代码分析

gffread gencode.v35.primary_assembly.annotation.gtf -g GRCh38.p13.genome.fa -w cds.fa --table @genename

在 cds.fa 数据中,我们就能得到 geneid 和 fasta 数据对应的文件,大概如下所示

head 5 cds.fa
>ENST00000456328.2	DDX11L1
GTTAACTTGCCGTCAGCCTTTTCTTTGACCTCTTCTTTCTGTTCATGTGTATTTGCTGTCTCTTAGCCCA
GACTTCCCGTGTCCTTTCCACCGGGCCTTTGAGAGGTCACAGGGTCTTGATGCTGTGGTCTTCATCTGCA
GGTGTCTGACTTCCAGCAACTGCTGGCCTGTGCCAGGGTGCAAGCTGAGCACTGGAGTGGAGTTTTCCTG
TGGAGAGGAGCCATGCCTAGAGTGGGATGGGCCATTGTTCATCTTCTGGCCCCTGTTGTCTGCATGTAAC
TTAATACCACAACCAGGCATAGGGGAAAGATTGGAGGAAAGATGAGTGAGAGCATCAACTTCTCTCACAA

接下来我们便可以进行读取,我这里读取的是每个转录本的最长的序列作为一个 geneid 的序列,这里我的 转录本id 和 geneid 之间的分隔符为 ' ',代码可以参考如下所示:

seq=protfasta.read_fasta("cds.fa", )
dic = {}
gene_names = []
for i in seq:
    k = i.split('\t')[-1]
    gene_names.append(k)
    try:
        if dic[k]:
            dic[k] = dic[k] if len(dic[k])>len(seq[i]) else seq[i]
            continue
    except:
        dic[k] = seq[i]
        continue

这样子便是一个得到一个字典,字典的 key 是 geneid,value 是其相对应的 fasta 序列,当我们想要再将其转为fasta文件时的时候,可以采取如下的代码:

f = open('gene_id.fasta','w')
for idx in zip(dic):
    i = idx[0]
    fa = dic[idx[0]]
    f.write('>'+i+'\n'+fa+'\n')
f.close()

查看文件可以看到

head gene_id.fasta
>gene0
ATGTCCCACCAACCGTCTATTATTCGAAAAACTCACCCTCTCCTATCATTAGGTAACAGTATATTAGTAGACCTCCCTTCTCCTGCTAACATCTCGGCCTGATGAAATTTTGGCTCATTATTAAGTTTATGCTTAATTCTACAAATTATTACTGGACTTATTCTTGCTATACACTACACCGCTAATACTGAACTAGCCTTCTCTTCAGTTATACACATTTGTCGTGACGTTAATAACGGATGACTTATACGAAACCTCCATGCTAATGGCGCCTCTATATTTTTTATCTGCATCTACGCTCATATCGGACGAGGAATTTATTATGGCTCCTATTTATATAAAGAAACATGAAACGTCGGAGTTATTTTATTTGCACTAACTGCAGCTACTGCCTTCGTAGGTTATGTTCTCCCATGGGGACAAATATCCTTTTGGGGGGCAACCGTTATCACAAATTTAATTTCAGCCATACCATATGTAGGAAATGATATTGTAGTATGATTATGGGGAGGCTTCTCAGTATCAAACGCCACTTTAACCCGATTCTTTACCTTCCATTTTATCTTACCATTCATTTTAGCAGCAATAACAATAATTCACATTATATTTCTTCACCAAACAGGATCTAGTAACCCTATAGGAATTAATTCTAATTTGGATAAGATTCAATTTCACCCGTATTTTTCTTTCAAAGATATTTTAGGTTTTGTTATTCTACTGGGCATTCTTTTCATAATTTCCCTTTTAGCCCCTAATGCACTAGGTGAACCAGACAACTTTATTTATGCTAATCCTCTTAGTACCCCTCCCCATATTAAACCAGAATGATACTTTCTATTTGCCTATGCCATTCTACGCTCTGTTCCTAATAAACTTGGAGGTGTTGTAGCTTTAGCAGCAGCTATCATAATCCTCCTAATTATCCCATTTACTCACACCTCCAAACAACGCGGAATACAATTTCGCCCACTCGCCCAAATTACATTTTGAATTTTAATTGCCGATCTAGCACTACTTACATGACTAGGGGGAGAGCCCGCTGAATATCCATTTATCTTAATAACACAAATTGCATCAACAGTCTACTTCATAATTTTTATTCTAGTTTTCCCAATTTTAGGATATTTAGAAAATAAAATACTATTAATATCAAAAAATACTGGTAAATTTAATTGAAAATTAGTTTACAGA