1 refseqID转换成geneID
下载gene2refseq.gz
wget -c ftp://ftp.ncbi.nih.gov/gene/DATA/gene2refseq.gz
2 提取相似物种的基因
gene2refseq.csv文件较大,有4G多,所以将blast比对过的物种中的提取出来,在进行ID转换。
首先查找物种的分类号
Glycine max 3847
Cicer arietinum 3827
Medicago truncatula 3880
vitis japonicus 29760
Homo sapiens 9606
Citrus sinensis 2711
Theobroma cacao 3641
Zea mays 4577-
grep命令(例)
grep ^9606 gene2refseq1 | grep ^38 gene2refseq1 | …
-
cut命令合并(应该不会写脚本)
cut gene2refseq* > gene.csv
(得到的两个文件:gene.csv和db_blast3.csv)
3 R语言-merge函数
merge函数的声明:
merge( x, y, by = intersect(names(x), names(y)),
by.x = by, by.y = by, all = FALSE, all.x = all, all.y = all,
sort = TRUE, suffixes = c(".x",".y"),
incomparables = NULL, ...)
#merge过程中两个文件匹配行一致,不然会报错。
DF1<-read.csv("lyr/rna-seq/blast/geneID/gene.csv") #读取gene.csv
DF2<-read.csv("lyr/rna-seq/blast/geneID/db_blast3.csv") #读取db_blast3.csv
dim(DF1) #看一下表格维度
dim(DF2)
merge(DF1,DF2,by ="RNA_nucleotide_accession.version", all.y = TRUE)
data1<-(merge(DF1,DF2,by ="RNA_nucleotide_accession.version", all.y = TRUE)) #将merge结果写道data中
write.csv(data1, file = "result_data.csv", quote = FALSE, row.names = FALSE) #输出文件为result_data.csv
接着看result_data.csv中基因的功能。