劉小澤寫于19.11.5
這次的推送要感謝小郭（幫我看到了一個(gè)非常容易忽視的錯(cuò)誤）和花花（幫我處理了一個(gè)簡(jiǎn)單但有用的腳本）
又讓我有了新的發(fā)現(xiàn)
問題是這樣滴
我的需求
我想下載SRA 數(shù)據(jù)，然后需要構(gòu)建一個(gè)配置文件（也就是一列是SRA ID）普气，一列是樣本名稱囤官。這樣做的目的是為了后面使用fastq-dump 進(jìn)行SRA轉(zhuǎn)fq文件的結(jié)果更易懂
# 比如原來(lái)的一個(gè)sra是SRR391032.sra扰魂。如果只使用默認(rèn)的fastq-dump參數(shù)，結(jié)果就是
SRR391032.1.fastq.gz
SRR391032.2.fastq.gz
# 但這樣的數(shù)據(jù)多了目代，我們就分不清哪個(gè)數(shù)據(jù)對(duì)應(yīng)哪個(gè)樣本屈梁，于是需要在轉(zhuǎn)換過程中就將樣本名對(duì)應(yīng)到fq文件上
# 例如構(gòu)建這樣一個(gè)config，第一列是SRR ID榛了，第二列是sample
# config文件就長(zhǎng)這樣(舉個(gè)例子)
SRR391032 WT-1
SRR391033 WT-2
SRR391034 TRT-1
SRR391035 TRT-2
現(xiàn)在找到GEO：https://www.ncbi.nlm.nih.gov//geo/query/acc.cgi?acc=GSE102339
進(jìn)來(lái)就能看到一些sample和GSM的對(duì)應(yīng)信息：
然后去獲得SRR和GSM的對(duì)應(yīng)信息：
https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP114984
打開之后是這樣（其中就包含了SRR和GSM的對(duì)應(yīng)信息）：
想法就是：將這兩個(gè)文件通過共同的GSM編號(hào)對(duì)應(yīng)起來(lái)
我的做法
step1：獲得GSM與樣本信息的對(duì)應(yīng)=》GSM-to-sample.txt
像??這樣：從Samples信息的底部直接復(fù)制到頂部在讶，然后粘貼到文本編輯器Sublime中，然后保存為GSM-to-sample.txt
然后讀到R中
step2：獲得GSM與SRR的對(duì)應(yīng)
也是將之前下載的SraRunTable.txt讀到R中
step3：準(zhǔn)備找二者的對(duì)應(yīng)關(guān)系
--------------------------- 分界線 -------------------------------
但就在這時(shí)霜大，我在mac上查看了step1的這個(gè)文件：GSM-to-sample.txt
less -SN GSM-to-sample.txt
看到是52行真朗，和網(wǎng)站上一致
強(qiáng)迫癥地想再次確認(rèn)一下行數(shù)，就出了問題
cat GSM-to-sample.txt| wc -l
# 51
這兩個(gè)竟然不一樣僧诚，差了一行遮婶！很奇怪蝗碎，想了很久，復(fù)制過來(lái)不應(yīng)該出錯(cuò)旗扑，而且less和讀到R中都是52行蹦骑。另外我單獨(dú)統(tǒng)計(jì)了第一列GSM號(hào)
cat GSM-to-sample.txt | cut -f1 | wc -l
# 52
小郭~復(fù)制粘貼少一行的解答
長(zhǎng)期從事Linux工作、具有豐富debug經(jīng)驗(yàn)臀防、孜孜不斷學(xué)習(xí)的小郭??一眼看出了問題：
”是不是因?yàn)樯倭艘粋€(gè)換行符“
真是一語(yǔ)點(diǎn)醒”夢(mèng)中人“（其實(shí)我不困??）
對(duì)啊眠菇，想起來(lái)我復(fù)制過來(lái)的最后一行好像真的沒有換行符，我重新vim編輯了一下這個(gè)文件袱衷，在最后一行敲了一下回車捎废，然后再刪掉最后的空行（這個(gè)看似沒有任何改變的操作，其實(shí)給最后一行加了一個(gè)換行符）
然后再次檢查：
cat GSM-to-sample.txt| wc -l
# 52
另外致燥，如果是linux環(huán)境登疗，還可以通過cat -A檢查最后時(shí)候具有換行符，會(huì)看到：
cat -A GSM-to-sample.txt
”為什么less命令查看就是52行呢嫌蚤？而wc統(tǒng)計(jì)就會(huì)不識(shí)別“
這個(gè)就可能是不同的命令識(shí)別機(jī)制了辐益，
wc -l統(tǒng)計(jì)的可能就是換行符（當(dāng)然我也沒有仔細(xì)去了解）
那么這個(gè)問題應(yīng)該如何避免呢？
可以直接去到Linux環(huán)境脱吱，先輸入一個(gè)cat >GSM-to-sample.txt
然后它會(huì)等待你的輸入智政，將復(fù)制的結(jié)果直接粘貼進(jìn)來(lái)，最后敲一下回車箱蝠，再ctrl + C 退出续捂，這時(shí)的文件也是完整的52行
總而言之，一定要確定最后一行有一個(gè)回車的操作宦搬，否則很有可能最后一行沒有換行符牙瓢，導(dǎo)致某些統(tǒng)計(jì)出錯(cuò)
花花~價(jià)值一元錢的腳本
現(xiàn)在有了兩個(gè)完整的文件，就要進(jìn)行SRR與樣本信息的對(duì)應(yīng)了
我當(dāng)時(shí)在思考上面的復(fù)制粘貼問題床三，為了省事，就花錢”雇“了花花幫我寫一個(gè)小腳本
腳本是這樣的：
rm(list = ls())
x1 =read.csv("~/Downloads/SRA-to-GSM.txt")
> x1[1:4,1:4]
Run Assay.Type AvgSpotLen BioProject
1 SRR6490544 ChIP-Seq 51 PRJNA397448
2 SRR5907429 ChIP-Seq 51 PRJNA397448
3 SRR5907430 ChIP-Seq 51 PRJNA397448
4 SRR5907431 ChIP-Seq 51 PRJNA397448
x2 = read.csv("~/Downloads/GSM-to-sample.txt",sep = "\t",header = F)
> head(x2,3)
GEO_Accession V2
1 GSM2734944 ChIP-Seq_Ez_WT_1
2 GSM2734945 ChIP-Seq_Ez_WT_2
3 GSM2734946 ChIP-Seq_Ez_WT_3
colnames(x2)[1] ="GEO_Accession"
x3 = merge(x2,x1,by = "GEO_Accession")
# 關(guān)于merge函數(shù)杨幼，花花的解釋是：”弱肉強(qiáng)食撇簿，誰(shuí)放在前面誰(shuí)是老大“，于是將x2放在前面就表示在x2的基礎(chǔ)上添加x1的內(nèi)容
> x3[1:4,1:4]
GEO_Accession V2 Run Assay.Type
1 GSM2734944 ChIP-Seq_Ez_WT_1 SRR5907429 ChIP-Seq
2 GSM2734945 ChIP-Seq_Ez_WT_2 SRR5907430 ChIP-Seq
3 GSM2734946 ChIP-Seq_Ez_WT_3 SRR5907431 ChIP-Seq
4 GSM2734947 ChIP-Seq_Pc_WT_1 SRR5907432 ChIP-Seq
config <- data.frame(srr=x3$Run,sample=x3$V2)
> head(config)
srr sample
1 SRR5907429 ChIP-Seq_Ez_WT_1
2 SRR5907430 ChIP-Seq_Ez_WT_2
3 SRR5907431 ChIP-Seq_Ez_WT_3
4 SRR5907432 ChIP-Seq_Pc_WT_1
5 SRR5907433 ChIP-Seq_Pc_WT_2
6 SRR5907434 ChIP-Seq_Ph_WT_1
沒想到這里又出現(xiàn)了一個(gè)BUG差购！
之前的GSM和樣本對(duì)應(yīng)關(guān)系是x1（經(jīng)上面的一頓操作四瘫，確實(shí)是52行）；GSM和SRR的對(duì)應(yīng)關(guān)系這里竟然只有51行（這個(gè)也沒錯(cuò)欲逃，的確是51行）找蜜。那么就說(shuō)明存在一個(gè)樣本不存在SRR
到底是哪個(gè)在搗亂？找出來(lái)稳析！
使用找差集的函數(shù)setdiff 洗做，依然符合之前花花說(shuō)的”弱肉強(qiáng)食“道理弓叛，誰(shuí)放在前面誰(shuí)是老大
這里的x2要多一行，所以它是老大
> setdiff(x2$GEO_Accession,x1$GEO_Accession)
[1] "GSM4030334"
在GEO官網(wǎng)找到這個(gè)GSM4030334：看一下https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM4030334
結(jié)果發(fā)現(xiàn)這么一句話：
因此诚纸，GEO上即使有GSM號(hào)撰筷，也不一定有SRR數(shù)據(jù)，所以做樣本信息與SRR對(duì)應(yīng)時(shí)畦徘，要小心了
歡迎關(guān)注我們的公眾號(hào)～_～
我們是兩個(gè)農(nóng)轉(zhuǎn)生信的小碩毕籽，打造生信星球，想讓它成為一個(gè)不拽術(shù)語(yǔ)井辆、通俗易懂的生信知識(shí)平臺(tái)关筒。需要幫助或提出意見請(qǐng)后臺(tái)留言或發(fā)送郵件到jieandze1314@gmail.com
