使用GTF和fasta生成每个geneid对应最长的fasta序列
在进行单细胞测序的时候,我们通常最后得到的是 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
All articles in this blog are licensed under CC BY-NC-SA 4.0 unless stating additionally.
Comment