搞定参考基因组，只需要五秒钟（序列相似性搜索工具

故事的开始是这样的~

看我如何用5个小时才解决了Jimmy老师5秒钟就帮我搞定的问题~
作为学徒的我这两天在跑Jimmy老师给的ATAC测试数据

ATAC测试数据

和chip-seq套路一样，第一步、第二步当然是：质控-->过滤，过滤完成之后使用bowtie2进行比对

这时我发现我不知道选哪个参考基因组比对，因为我根本不知道这个测试数据是什么物种？

以前跑流程的时候用的是人家文章里的数据，测序数据是人是鼠、比对哪个参考基因组，文章里写的很清楚，所以这个问题以前根本就没在意过。在江湖上混，曾经没学会的东西迟早要还的，好吧，现在问题来了：

在没有文章只有测序数据的情况下我们怎么知道序列是人是鼠？

询问健明老师后，老师给了我简单、明确、直击要点的回答：

其实，刚看到blat这四个英文字母的时候我发现我脑子里没有关于它的内存，经过我苦苦搜索，终于找到了UCSC上的网页工具BLAT，我以为我的问题马上就得到解决了，然而……

知道了UCSC上有这么个神奇的网页工具后，我就马不停蹄的去谷歌上搜索

UCSC

看到一个顺眼的网页链接就直接点进去了。

胜利在望的喜悦冲昏了我的本来就不灵光的头脑，蒙蔽了我本来就近视的双眼zea mays

那两大字被我无情忽视，直接把序列放在框里点submit提交了

结果就一直是找不到，在我复制了二三十条序列之后还一直出现这个结果，我根本就没意识到是我进错网页啦。

其实BLAT这个网页小工具在健明老师的B站视频讲过的，如果我早些看过就不会掉这个坑里并在坑里挣扎5个小时了……我可以直接去跳下一个坑啦，前面那么多坑等着我往里跳呢，干嘛在这个坑里浪费时间~

下面我们看看如何得到fastq文件里的序列信息，看下我的测试数据：

$ ls
C1_R1.fastq.gz C1_R2.fastq.gz

是两个fastq压缩文件，再看下fastq文件的格式：四行代表一条reads信息，第二行就是我们想要的序列

fastq文件格式用zcat查看压缩文件，直接把第二行复制到BLAT里那个放序列的框里就行

zcat C1_R1.fastq.gz | less -S

$ zcat C1_R1.fastq.gz | less -S
@E00515:437:HTN33CCXY:3:1101:12276:1731 1:N:0:TCCTTCCA
ACATAAGCTGAAAGACATCTAGACCTGGATCTGTTGTGCTCACCGTTGTCTTGATCATGTTTCAGGATTGATGCTAAATAAATATCTTGAGTGCATGGATGATT
+
AAF-AJJJJAJJFJAJ<FJJJ<-<-<FJ-F<-J<A-<F-<JJJJJJJJ--J-F--F--F-F-<----<<7<A-AF-FFA--AF-AJ<FJFA-A777---7-<77
@E00515:437:HTN33CCXY:3:1101:12560:1731 1:N:0:TCCTTCCA
GCCACACAGATGCCAGGGCAGATTCAGGAAGGGTCGTGGCCTCATGGGATCGATCAGTGACGCATCTTTGAGGCCACTAGCAGCCCGTGTCTCTTATACGCATC
+
A-AA<FJJJFJJJJFJJJJJJJJF-FJJ--F-JAJ-7JFJJJJJFJJJA-FAJ-F<-7<-F7AJ----FFFJJJJAJ<AJF7F-7J7-A-FF7-<7--7-<7--
@E00515:437:HTN33CCXY:3:1101:16396:1731 1:N:0:TCCTTCCA
AGTTTACACCCCTTGTGTGACTTTGCTCCTAAAGGGGAGCGGGGGAGAGCAAGATAGTAGAGCCGCTATCGCAGGGGGTGTACCTGTCTCTTATACAGATCTCG
+
AAF-AFFJJAJAFFJFJFJJJ<FF7JFF--7-FJJJ<7FAJJJJJ<JFF-A7J7<J--7<<--7--7-<-FFFJ-AFJ-7-FJF7JF---7---AA-----7--
@E00515:437:HTN33CCXY:3:1101:18588:1731 1:N:0:TCCTTCCA
CAAATAAAGAATAGGCGATAGTAGTTGAAACGTGGGGCAATAGATATAGTACCTCAAGGCATAGATGAAAAATTATAACCAAGCATAATATAGCAAGGTCTAAC
+
-AAFFJJJJFJJ-FJJJAJJJ<J77FJJ7JF-FFJF-AJJJJJJJJFJJJJAJAJA-FJ<F-JJ--<AFAJJJJAFJJ-JFJJAJJ7JJJJF<AFAA-7-A<7A
@E00515:437:HTN33CCXY:3:1101:21775:1731 1:N:0:TCCTTCCA
GCTCAAGAATATGCTGTCCTAGTGGGTGACAGAGGTCAGGCAGGAATTAAGCAGGTGTTGCTATGTAAAAGGCTGAGCGTACAGTAGGTTTGCACCAGGGATAG

如果我看其他三行碍事，影响了我复制粘贴，只想要第二行序列怎么办？好办
zcat C1_R1.fastq.gz | paste - - - - | cut -f2 |less -S

$ zcat C1_R1.fastq.gz | paste - - - - | cut -f2 |less -S
ACATAAGCTGAAAGACATCTAGACCTGGATCTGTTGTGCTCACCGTTGTCTTGATCATGTTTCAGGATTGATGCTAAATAAATATCTTGAGTGCATGGATGATT
GCCACACAGATGCCAGGGCAGATTCAGGAAGGGTCGTGGCCTCATGGGATCGATCAGTGACGCATCTTTGAGGCCACTAGCAGCCCGTGTCTCTTATACGCATC
AGTTTACACCCCTTGTGTGACTTTGCTCCTAAAGGGGAGCGGGGGAGAGCAAGATAGTAGAGCCGCTATCGCAGGGGGTGTACCTGTCTCTTATACAGATCTCG
CAAATAAAGAATAGGCGATAGTAGTTGAAACGTGGGGCAATAGATATAGTACCTCAAGGCATAGATGAAAAATTATAACCAAGCATAATATAGCAAGGTCTAAC
GCTCAAGAATATGCTGTCCTAGTGGGTGACAGAGGTCAGGCAGGAATTAAGCAGGTGTTGCTATGTAAAAGGCTGAGCGTACAGTAGGTTTGCACCAGGGATAG
TTTTTACCCTATGCAATTTCTAATGGTTGGAGAAACCAAGGCTGCCAGTTCCATC
paste - - - -`的意思是，把四行剪切粘贴成一行，把四行粘成一行后，行与行之间默认分隔符为 Tab . 然后 `cut -f2`取第二个字段便是序列信息。
这代码写的太low了，能不能有一个高级一点的写法来衬托一下我已经从小白变成菜鸟了？有
`zcat C1_R1.fastq.gz | awk '{if(NR%4==2)print}' | less -S
$ zcat C1_R1.fastq.gz | awk '{if(NR%4==2)print}' | less -S
ACATAAGCTGAAAGACATCTAGACCTGGATCTGTTGTGCTCACCGTTGTCTTGATCATGTTTCAGGATTGATGCTAAATAAATATCTTGAGTGCATGGATGATT
GCCACACAGATGCCAGGGCAGATTCAGGAAGGGTCGTGGCCTCATGGGATCGATCAGTGACGCATCTTTGAGGCCACTAGCAGCCCGTGTCTCTTATACGCATC
AGTTTACACCCCTTGTGTGACTTTGCTCCTAAAGGGGAGCGGGGGAGAGCAAGATAGTAGAGCCGCTATCGCAGGGGGTGTACCTGTCTCTTATACAGATCTCG
CAAATAAAGAATAGGCGATAGTAGTTGAAACGTGGGGCAATAGATATAGTACCTCAAGGCATAGATGAAAAATTATAACCAAGCATAATATAGCAAGGTCTAAC

NR是awk里面的一个内置变量，代表文件行号。
NR%4==2就是取第二行的意思
你骗人，这个代码也不好看，还有没有更好的？好吧，我承认，我没有啦！因为我还是一个小白还没有变成菜鸟。想要更好的就看下Jimmy老师在B站上的Linux视频吧

如果你觉光看视频还是不够，那就来跟我一起当学徒吧，享受VIP一对一指导。

下面是正文

今天我们介绍一款网页版序列相似性搜索工具—UCSC BLAT，这款小工具能方便快捷的告诉我们目标序列是人的还是小鼠的。

操作起来特别简单：

1. 打开BLAT主页

2. 选择合适参数进行序列搜索

3. 查看序列搜索结果

干货部分到此结束！

只讲干货不讲实际应用以及我是怎么掉坑里的，那不是我的风格~~~

搞定参考基因组，只需要五秒钟（序列相似性搜索工具—UCSC BLAT）

相关推荐