sample等价是什么错误_一个复制粘贴引发的有趣小错误及思考

今天是生信星球陪你的第481天

大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~

就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~

这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！

豆豆写于19.11.5
这次的推送要感谢小郭(帮我看到了一个非常容易忽视的错误)和花花(帮我处理了一个简单但有用的脚本)
又让我有了新的发现

问题是这样滴

我的需求

我想下载SRA 数据，然后需要构建一个配置文件(也就是一列是SRA ID)，一列是样本名称。这样做的目的是为了后面使用fastq-dump 进行SRA转fq文件的结果更易懂

# 比如原来的一个sra是SRR391032.sra。如果只使用默认的fastq-dump参数，结果就是SRR391032.1.fastq.gzSRR391032.2.fastq.gz# 但这样的数据多了，我们就分不清哪个数据对应哪个样本，于是需要在转换过程中就将样本名对应到fq文件上# 例如构建这样一个config，第一列是SRR ID，第二列是sample# config文件就长这样(举个例子)SRR391032 WT-1SRR391033 WT-2SRR391034 TRT-1SRR391035 TRT-2

现在找到GEO：https://www.ncbi.nlm.nih.gov//geo/query/acc.cgi?acc=GSE102339

进来就能看到一些sample和GSM的对应信息：

然后去获得SRR和GSM的对应信息：

https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP114984

打开之后是这样(其中就包含了SRR和GSM的对应信息)：

想法就是：将这两个文件通过共同的GSM编号对应起来

我的做法

step1：获得GSM与样本信息的对应=》`GSM-to-sample.txt`

像?这样：从Samples信息的底部直接复制到顶部，然后粘贴到文本编辑器Sublime中，然后保存为GSM-to-sample.txt

然后读到R中

step2：获得GSM与SRR的对应

也是将之前下载的SraRunTable.txt读到R中

step3：准备找二者的对应关系

--------------------------- 分界线 -------------------------------

但就在这时，我在mac上查看了step1的这个文件：GSM-to-sample.txt

less -SN GSM-to-sample.txt

看到是52行，和网站上一致

强迫症地想再次确认一下行数，就出了问题

cat GSM-to-sample.txt| wc -l# 51

这两个竟然不一样，差了一行！很奇怪，想了很久，复制过来不应该出错，而且less和读到R中都是52行。另外我单独统计了第一列GSM号

cat GSM-to-sample.txt | cut -f1 | wc -l# 52

小郭~复制粘贴少一行的解答

长期从事Linux工作、具有丰富debug经验、孜孜不断学习的小郭?一眼看出了问题：

”是不是因为少了一个换行符“

真是一语点醒”梦中人“(其实我不困?)

对啊，想起来我复制过来的最后一行好像真的没有换行符，我重新vim编辑了一下这个文件，在最后一行敲了一下回车，然后再删掉最后的空行(这个看似没有任何改变的操作，其实给最后一行加了一个换行符)

然后再次检查：

cat GSM-to-sample.txt| wc -l# 52

另外，如果是linux环境，还可以通过cat -A检查最后时候具有换行符，会看到：

cat -A GSM-to-sample.txt

”为什么less命令查看就是52行呢？而wc统计就会不识别“

这个就可能是不同的命令识别机制了，wc -l统计的可能就是换行符(当然我也没有仔细去了解)

那么这个问题应该如何避免呢？

可以直接去到Linux环境，先输入一个cat >GSM-to-sample.txt

然后它会等待你的输入，将复制的结果直接粘贴进来，最后敲一下回车，再ctrl + C 退出，这时的文件也是完整的52行

总而言之，一定要确定最后一行有一个回车的操作，否则很有可能最后一行没有换行符，导致某些统计出错

花花~价值一元钱的脚本

现在有了两个完整的文件，就要进行SRR与样本信息的对应了

我当时在思考上面的复制粘贴问题，为了省事，就花钱”雇“了花花帮我写一个小脚本

脚本是这样的：

rm(list = ls())x1 =read.csv("~/Downloads/SRA-to-GSM.txt") > x1[1:4,1:4]         Run Assay.Type AvgSpotLen  BioProject1 SRR6490544   ChIP-Seq         51 PRJNA3974482 SRR5907429   ChIP-Seq         51 PRJNA3974483 SRR5907430   ChIP-Seq         51 PRJNA3974484 SRR5907431   ChIP-Seq         51 PRJNA397448

x2 = read.csv("~/Downloads/GSM-to-sample.txt",sep = "\t",header = F)> head(x2,3)  GEO_Accession               V21    GSM2734944 ChIP-Seq_Ez_WT_12    GSM2734945 ChIP-Seq_Ez_WT_23    GSM2734946 ChIP-Seq_Ez_WT_3

colnames(x2)[1] ="GEO_Accession"x3 = merge(x2,x1,by = "GEO_Accession")# 关于merge函数，花花的解释是：”弱肉强食，谁放在前面谁是老大“，于是将x2放在前面就表示在x2的基础上添加x1的内容> x3[1:4,1:4]  GEO_Accession               V2        Run Assay.Type1    GSM2734944 ChIP-Seq_Ez_WT_1 SRR5907429   ChIP-Seq2    GSM2734945 ChIP-Seq_Ez_WT_2 SRR5907430   ChIP-Seq3    GSM2734946 ChIP-Seq_Ez_WT_3 SRR5907431   ChIP-Seq4    GSM2734947 ChIP-Seq_Pc_WT_1 SRR5907432   ChIP-Seq

config > head(config)         srr           sample1 SRR5907429 ChIP-Seq_Ez_WT_12 SRR5907430 ChIP-Seq_Ez_WT_23 SRR5907431 ChIP-Seq_Ez_WT_34 SRR5907432 ChIP-Seq_Pc_WT_15 SRR5907433 ChIP-Seq_Pc_WT_26 SRR5907434 ChIP-Seq_Ph_WT_1

没想到这里又出现了一个BUG！

之前的GSM和样本对应关系是x1(经上面的一顿操作，确实是52行)；GSM和SRR的对应关系这里竟然只有51行(这个也没错，的确是51行)。那么就说明存在一个样本不存在SRR

到底是哪个在捣乱？找出来！

使用找差集的函数setdiff ，依然符合之前花花说的”弱肉强食“道理，谁放在前面谁是老大

这里的x2要多一行，所以它是老大

> setdiff(x2$GEO_Accession,x1$GEO_Accession)[1] "GSM4030334"

在GEO官网找到这个GSM4030334：看一下https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM4030334

结果发现这么一句话：

因此，GEO上即使有GSM号，也不一定有SRR数据，所以做样本信息与SRR对应时，要小心了

初学生信，很荣幸带你迈出第一步

?生信星球 ?~ 一个不拽术语、通俗易懂的生信知识平台