故事的开始是这样的~
看我如何用5个小时才解决了Jimmy老师5秒钟就帮我搞定的问题~ 作为学徒的我这两天在跑Jimmy老师给的ATAC测试数据
ATAC测试数据
和chip-seq套路一样,第一步、第二步当然是:质控-->过滤,过滤完成之后使用bowtie2进行比对
这时我发现我不知道选哪个参考基因组比对,因为我根本不知道这个测试数据是什么物种?
以前跑流程的时候用的是人家文章里的数据,测序数据是人是鼠、比对哪个参考基因组,文章里写的很清楚,所以这个问题以前根本就没在意过。在江湖上混,曾经没学会的东西迟早要还的,好吧,现在问题来了:
其实,刚看到blat这四个英文字母的时候我发现我脑子里没有关于它的内存,经过我苦苦搜索,终于找到了UCSC上的网页工具BLAT,我以为我的问题马上就得到解决了,然而……
知道了UCSC上有这么个神奇的网页工具后,我就马不停蹄的去谷歌上搜索
代码语言:javascript复制UCSC
看到一个顺眼的网页链接就直接点进去了。
胜利在望的喜悦冲昏了我的本来就不灵光的头脑,蒙蔽了我本来就近视的双眼zea mays
那两大字被我无情忽视,直接把序列放在框里点submit提交了
结果就一直是找不到,在我复制了二三十条序列之后还一直出现这个结果,我根本就没意识到是我进错网页啦。
其实BLAT这个网页小工具在健明老师的B站视频讲过的,如果我早些看过就不会掉这个坑里并在坑里挣扎5个小时了……我可以直接去跳下一个坑啦,前面那么多坑等着我往里跳呢,干嘛在这个坑里浪费时间~
下面我们看看如何得到fastq
文件里的序列信息,看下我的测试数据:
$ 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
$ 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. 查看序列搜索结果
干货部分到此结束!
只讲干货不讲实际应用以及我是怎么掉坑里的,那不是我的风格~~~