EST:表达序列标签,expressed sequence tags 。
顾名思义,很好理解,就是表达出来的序列,即从基因组DNA上表达出来的RNA,但是我们没法测序RNA,所以我们最终测的是表达序列的cDNA片段。
“标签”:就是指这些序列可能比较短,但是可以用来标定一个物种。
常见下载方式有两种:
1. NCBI Web下载
打开,搜索你要的物种,比如 Camellia ,可以看到结果EST (50287)。
Web下载几个还行,想要批量下载就有点费力了,ncbi反爬虫,也不好爬。
2. NCBI ftp下载
直接wget就可以批量下载了。
for one in `seq 1 81`doecho $onewget ftp://ftp.ncbi.nih.gov/repository/dbEST/dbEST.reports.000000.${one}.gzdone
下载后的文件格式是:
IDENTIFIERSdbEST Id: 5EST name: EST00006GenBank Acc: M61958GDB Dsegment: D0S2525ECLONE INFOClone Id: HHCSB86Source: ATCCId in host: 77063DNA type: cDNAPRIMERSSequencing: M13 ForwardPolyA Tail: UnknownSEQUENCE TGCACAACCAAGTTTTGTGACTACGGGAAGGCTCCCGGGGCAGAGGAGTACGCTCAACAA GATGTGTTAAAGAAATCTTACTCCAAGGCCTTCACGCTGACCATCTCTGCCCTCTTTGTG ACACCCAAGACGACTGGGGCCCNGGTGGAGTTAAGCGAGCAGCAACTNCAGTTGTNGCCG AGTGATGTGGACAAGCTGTCACCCACTGACAEntry Created: May 26 1992Last Updated: Dec 18 2012PUTATIVE ID Assigned by submitter 2',3'-cyclic nucleotide phoshodiesteraseLIBRARYId: LIBEST_000004Lib Name: LIBEST_000004 Hippocampus, Stratagene (cat. #936205)Organism: Homo sapiensVector: lambdaZAP-IIDescription: Female, 2 years; oligo-dT + random primed cDNA synthesis;
信息是挺全面的,自己想要哪个物种就只能自己提取了。
提取成FASTA的脚本我就不贴了(效率很重要,因为文件很大)。
最后我还是自己写了个脚本,biopython实在是太慢了。
import gzipinf = gzip.open("dbEST.reports.000000.49.gz","rb")raw_id = ""seq = ""for line in inf: if line.stratswith("GenBank Acc"): id = line.split(":")[1].strip() if line.stratswith("SEQUENCE"): seq = "" while True: rline = inf.readline() seq+=rline.strip() if not line.stratswith(" "): break if line.stratswith("Organism"): organism = line.split(":")[1].strip() if organism.startswith("Camellia"): print(">"+id+" "+organism, seq, sep="\n")
我的脚本可以用,但是不一定很快。
我用awk试了很久,没有成功。
2018年3月16日