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

2019-05-14 18:49:04 浏览数 (1)

故事的开始是这样的~

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

ATAC测试数据

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

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

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

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

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

代码语言:javascript复制
UCSC

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

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

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

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

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


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

代码语言:javascript复制
$ ls
C1_R1.fastq.gz  C1_R2.fastq.gz

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

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

代码语言:javascript复制
zcat C1_R1.fastq.gz | less -S
代码语言:javascript复制
$ 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

代码语言:javascript复制
$ 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. 查看序列搜索结果

干货部分到此结束!

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

0 人点赞