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

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

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

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

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

1
2
3
4
5
6
7
head 5 cds.fa
>ENST00000456328.2	DDX11L1
GTTAACTTGCCGTCAGCCTTTTCTTTGACCTCTTCTTTCTGTTCATGTGTATTTGCTGTCTCTTAGCCCA
GACTTCCCGTGTCCTTTCCACCGGGCCTTTGAGAGGTCACAGGGTCTTGATGCTGTGGTCTTCATCTGCA
GGTGTCTGACTTCCAGCAACTGCTGGCCTGTGCCAGGGTGCAAGCTGAGCACTGGAGTGGAGTTTTCCTG
TGGAGAGGAGCCATGCCTAGAGTGGGATGGGCCATTGTTCATCTTCTGGCCCCTGTTGTCTGCATGTAAC
TTAATACCACAACCAGGCATAGGGGAAAGATTGGAGGAAAGATGAGTGAGAGCATCAACTTCTCTCACAA

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

1
2
3
4
5
6
7
8
9
10
11
12
13
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文件时的时候，可以采取如下的代码：

1
2
3
4
5
6
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()

查看文件可以看到

1
2
3
head gene_id.fasta
>gene0
ATGTCCCACCAACCGTCTATTATTCGAAAAACTCACCCTCTCCTATCATTAGGTAACAGTATATTAGTAGACCTCCCTTCTCCTGCTAACATCTCGGCCTGATGAAATTTTGGCTCATTATTAAGTTTATGCTTAATTCTACAAATTATTACTGGACTTATTCTTGCTATACACTACACCGCTAATACTGAACTAGCCTTCTCTTCAGTTATACACATTTGTCGTGACGTTAATAACGGATGACTTATACGAAACCTCCATGCTAATGGCGCCTCTATATTTTTTATCTGCATCTACGCTCATATCGGACGAGGAATTTATTATGGCTCCTATTTATATAAAGAAACATGAAACGTCGGAGTTATTTTATTTGCACTAACTGCAGCTACTGCCTTCGTAGGTTATGTTCTCCCATGGGGACAAATATCCTTTTGGGGGGCAACCGTTATCACAAATTTAATTTCAGCCATACCATATGTAGGAAATGATATTGTAGTATGATTATGGGGAGGCTTCTCAGTATCAAACGCCACTTTAACCCGATTCTTTACCTTCCATTTTATCTTACCATTCATTTTAGCAGCAATAACAATAATTCACATTATATTTCTTCACCAAACAGGATCTAGTAACCCTATAGGAATTAATTCTAATTTGGATAAGATTCAATTTCACCCGTATTTTTCTTTCAAAGATATTTTAGGTTTTGTTATTCTACTGGGCATTCTTTTCATAATTTCCCTTTTAGCCCCTAATGCACTAGGTGAACCAGACAACTTTATTTATGCTAATCCTCTTAGTACCCCTCCCCATATTAAACCAGAATGATACTTTCTATTTGCCTATGCCATTCTACGCTCTGTTCCTAATAAACTTGGAGGTGTTGTAGCTTTAGCAGCAGCTATCATAATCCTCCTAATTATCCCATTTACTCACACCTCCAAACAACGCGGAATACAATTTCGCCCACTCGCCCAAATTACATTTTGAATTTTAATTGCCGATCTAGCACTACTTACATGACTAGGGGGAGAGCCCGCTGAATATCCATTTATCTTAATAACACAAATTGCATCAACAGTCTACTTCATAATTTTTATTCTAGTTTTCCCAATTTTAGGATATTTAGAAAATAAAATACTATTAATATCAAAAAATACTGGTAAATTTAATTGAAAATTAGTTTACAGA