一个朋友的文章被审稿人质疑数据真实性，如何凭实力Argue？（明明这个数据集的平台注释文件没有lncRNA吖，你的数据集怎么会有的，我对你的数据真实性持严重怀疑态度）
数据集 GSE65858
平台文件 GPL10558
其实就是对原来的平台注释文件做个重注释，重新比对到转录本上。
准备工作
1.转录本信息文件下载
进入genecode数据库，选择“人”，最后选择最新的比对文件：gencode.v38.transcripts.fa.gz
demo1
2.GEO平台文件下载 (有点大，39G)
demo2
3.比对软件下载
seqmap是用于短序列比对特别快的工具。但是它出来的结果没有blast和blat多，如果要对芯片的序列进行重注释，是很好的一个工具.
搜索引擎输入seqmap，安装适合自己机器的版本就行了
这里因为平台信息文件太大了，所以，我在Linux下完成这些操作。
如果是window系统，可以直接来这个seqmap下载。
如果是Linux系统：
conda install -c bioconda seqmap
conda install -c bioconda/label/cf201901 seqmap
seqmap常规参数
软件的基本格式为：seqmap [options]
1.输入格式中参考基因组和比对的基因组必须是fa格式
2.num_mismatch代表比对的时候不匹配的个数
3.输出文件的格式分为两种。其中默认的是：Eland格式。另外一种是我们可以看得比较清楚的。用来显示所有匹配结果的格式：/output_all_matches
4.处理平台文件变成fasta
最困难的地方是fasta文件获取，这个只能从平台文件来转换。（这里其实也可以用平台注释文件代替，如果里面包括了sequence序列信息的话）
gpl <- data.table::fread("GPL10558_family.soft",skip = "ID",data.table = F)
## 可以使用平台注释文件代替（包含了sequence）
gpl <- data.table::fread("GPL10558-50081.txt",skip = "ID",data.table = F)
## 现在就选择两列，一列是探针名称，一列是序列
## 选取想要的两列，一列是ID， 一列是序列
gpl <- gpl[,c(1,22)]
## 如果选用了平台注释文件
gpl <- gpl[,c(1,19)]
## 过滤掉没有序列的行
library(dplyr)
gpl <- gpl %>%
filter(nchar(SEQUENCE)!=0)
## 保存为fasta格式文件,十分巧妙！
gp <- paste0('>',gpl$ID,'\n', gpl$SEQUENCE)
## 写成文本文件的时候就可以看出fasta的格式了
write.table(gp,'GPL.fasta', quote = F, row.names = F, col.names = F)
此时，我们要的三个文件都齐了
1.序列比对文件
2.参考基因组
3.seqmap软件
5.批量比对
./seqmap-1.0.12-linux 0 ./GPL.fasta ./gencode.v38.transcripts.fa seqmap_results.txt /output_all_matches
这个命令由6部分组成
./seqmap-1.0.12-windows.exe是软件名称，
0表示匹配容错率为0
./GPL.fasta是平台fasta文件
./gencode.v30.transcripts.fa是参考基因组
seqmap_results.txt是生成文件的名字
//output_all_matchest输出所有匹配结果，其他系统一个斜杠即可。
比对
6.提取比对结果
## 读入数据
probe2ID <- data.table::fread("seqmap_results.txt",data.table = F)
## 重要的结果在第一列
library(tidyr)
library(dplyr)
probe2ID <- probe2ID %>%
select(probe_id,trans_id) %>%
separate(trans_id,into = c("Ensembl",
"drop1","drop2","drop3",
"trans_Symble","gene_Symble","drop4","trans_biotype"),sep = "\\|") %>%
select(probe_id,Ensembl,trans_Symble,gene_Symble,trans_biotype)
最终结果如下图：
比对结果
网友评论