一个引号引发的血案

Posted on by

安排学徒探索了一下表达量芯片的不同探针平台信息，然后学徒给我反馈了一个在他看来有意思的bug，就是在读取一个txt文件的时候会出现读不完整的情况 ：

k = read.table('./GPL570-55999.txt',
 header = T,
 sep = '\t',
 fill = T,
 #quote = '',
 check.names = F,
 skip = 16)

Warning message:
In scan(file = file, what = what, sep = sep, quote = quote, dec = dec, :
 EOF within quoted string

Warning 信息也得注意

这个时候函数给出来的是Warning message，并不是大家常见的error，所以有些人可能会忽略它。

其实解决方式很简单，需要仔细调整参数，比如在上面的read.table 函数里面添加了 quote = ‘’ 的参数，或者使用高级的R包 ：

k2 = rio::import('./GPL570-55999.txt',skip = 16)
k1 = data.table::fread('./GPL570-55999.txt',skip = 16)

上面的失败读取的问题在于 最后一个探针 “AFFX-r2-Bs-dap-M_at” ：

> dim(k);dim(k1);dim(k2)
[1] 28059 16
[1] 54675 16
[1] 54675 16

> head(k[,1]);tail(k[,1])
[1] "1007_s_at" "1053_at" "117_at" "121_at" "1255_g_at"
[6] "1294_at" 
[1] "AFFX-M27830_M_at" "AFFX-PheX-3_at" "AFFX-PheX-5_at" 
[4] "AFFX-PheX-M_at" "AFFX-r2-Bs-dap-3_at" "AFFX-r2-Bs-dap-M_at"

所以需要去看看这个探针对应的注释信息，这个时候不用R代码了，就一个简简单单shell命令即可 ：

grep -n 'AFFX-r2-Bs-dap-M_at' GPL570-55999.txt 
54668:AFFX-r2-Bs-dap-M_at --Control Homo sapiens Oct 6, 2014 Control sequence Affymetrix Proprietary Database B. subtilis /GEN=dapB, jojF /DB_XREF=gb:L38424.1 /NOTE=SIF corresponding to nucleotides 2055-2578 of gb:L38424.1 /DEF=Bacillus subtilis dihydropicolinate reductase (jojE) gene, complete cds; poly(A) polymerase (jojI) gene, complete cds; biotin acetyl-CoA-carboxylase ligase (birA) gene, complete cds; jojC, jojD, jojF, jojG, jojH genes, complete cds's. AFFX-r2-Bs-dap-M
-ne

其实我本来是准备让搜索如何输出指定的第 28059 行 附近的信息，还使用人工智能大模型帮我搜索了方法 （在Linux的shell里面如何针对txt文件输出指定的某一行）：

在Linux的shell中，可以使用几种不同的方法来输出文本文件中的指定行。以下是一些常用的命令 ：

  1. 使用sed命令：

    • sed 'Nq' filename.txt：输出文件filename.txt的第N行。

  2. 使用awk命令：

    • awk 'NR==N {print}' filename.txt：输出文件filename.txt的第N行，其中NR是awk内置的记录数变量。

  3. 使用head和tail命令（当你知道行号在文件的前几行时）：

    • head -n N filename.txt | tail -n 1：输出文件filename.txt的第N行，适用于行号较小的情况。

  4. 使用grep命令（如果行包含特定的文本）：

    • grep 'pattern' filename.txt：输出包含特定文本pattern的所有行。使用-n选项可以显示行号。

  5. 使用less或more命令（交互式查看）：

    • 可以先使用less或more命令查看文件内容，然后使用上下箭头键跳转到指定的行。

  6. 使用nl命令（带行号查看）：

    • nl filename.txt：显示文件filename.txt的所有行，并带上行号。然后你可以手动查找到指定的行。

  7. 使用cat与-n选项：

    • cat -n filename.txt：类似于nl，cat -n命令也会输出文件的每一行，并在行首显示行号。

  8. 使用tac命令（从文件末尾开始计数）：

    • tac filename.txt | grep '^[N]'：这将从文件末尾开始计数，显示倒数第N行。注意，^在正则表达式中表示行的开始，[N]是十进制数。

选择哪种方法取决于你的具体需求和文件的大小。对于大型文件，建议使用awk或sed，因为它们在处理大文件时更高效。对于较小的文件，或者当你只需要查看文件的一部分时，head、tail或grep可能更简单快捷。

首先，这个大模型给出来的第一个方法就是错误的，然后呢我自己提出来的输出指定的第 28059 行 附近的信息这个思路是错误的，因为read.table 函数读取txt后的数据框的大小并不是txt文件本身出问题的地方。因为我们选择 在上面的read.table 函数里面添加了 quote = ‘’ 的参数，其实就明白了是因为txt文件里面有单引号，这个才是引发问题的原因，文本文件里面的有 1333 个地方都有这个单引号 ：

grep "'" GPL570-55999.txt |wc
 1333

有些时候连Warning都没有但也是错误的

比如另外一个GPL的txt文件读取：

k2 = rio::import('./GPL13158-5065.txt',skip = 16)
k1 = data.table::fread('./GPL13158-5065.txt',skip = 16)
k = read.table('./GPL13158-5065.txt',
 header = T,
 sep = '\t',
 fill = T,
 check.names = F,
 skip = 16)

完全是顺利的运行了：

> dim(k);dim(k1);dim(k2)
[1] 26374 16
[1] 54715 16
[1] 54715 16
> head(k[,1]);tail(k[,1])
[1] "1007_PM_s_at" "1053_PM_at" "117_PM_at" "121_PM_at" 
[5] "1255_PM_g_at" "1294_PM_at" 
[1] "AFFX-ThrX-3_at" "AFFX-ThrX-5_at" "AFFX-ThrX-M_at" 
[4] "AFFX-TrpnX-3_at" "AFFX-TrpnX-5_at" "AFFX-TrpnX-M_at"
>

高级函数也会报错

前面提到了用高级的R包 ：

k2 = rio::import('./GPL570-55999.txt',skip = 16)
k1 = data.table::fread('./GPL570-55999.txt',skip = 16)

但是它们也不是万能的哦，之前就遇到了一个单细胞转录组表达量矩阵文件，是txt或者csv格式的， 使用fread就只能读取一半的基因或者细胞。

​

Comments are closed.