有这个需求，是因为我们经常对某些细胞系进行有针对性的设计变异，比如BAF155的R1064K呀，H3F3A的K27呀，那我我们拿到高通量测序数据的时候，就肯定希望可以快速的看看这个基因是否被突变成功了。现在比对几乎不耗费什么时间了，但是得到的sam要sort的时候还是蛮耗费时间的。假设，我们已经得到了所有样本的sort好的bam文件，想看看自己设计的基因突变是否成功了，可以有针对性的只call 某个基因的突变！
十一
14
Monthly Archives: 11月 2016
十一
12
仔细探究picard的MarkDuplicates 是如何行使去除PCR重复reads功能的
本帖紧跟前面的仔细探究samtools的rmdup是如何行使去除PCR重复reads功能的
同样的我们也是分单端和双端测序来看结果，并且比较两个工具的区别！
首先对于那个单端数据，samtools给出的结果是：[bam_rmdupse_core] 25 / 53 = 0.4717 in library Continue reading
十一
12
仔细探究samtools的rmdup是如何行使去除PCR重复reads功能的
在做这个去除PCR重复reads时候必须要明白为什么要做这个呢？WGS？WES？RNA-SEQ?CHIP-SEQ？都需要吗？随机打断测序才需要？特异性捕获不需要？
搞明白了，我们就开始做，首先拿一个小的单端测序数据比对结果来做测试！
samtools rmdup -s tmp.sorted.bam tmp.rmdup.bam
[bam_rmdupse_core] 25 / 53 = 0.4717 in library
我们的测试数据里面有53条records根据软件算出了25条reads都是PCR的duplicate，所以去除了！
十一
11
数据库批量注释不可盲目-annovar数据库错误
我对H3F3A这个基因做了两个突变的cellline，分别是G34V和K27M，现在知道这个基因在hg38上面的坐标是：
Genomic Location for H3F3A Gene
Chromosome: 1
Start:226,061,851 bp from pter End:226,072,002 bp from pter
Size:10,152 bases Orientation:Plus strand
然后我用samtools结合bcftools把该基因区域的snp位点call出来：
samtools mpileup -r chr1:226061851-226072001 -t "DP4" -ugf ~/reference/genome/hg38/hg38.fa *sorted.bam | bcftools call -vmO z -o H3F3A.vcf.gz
十一
10
一个基因在同一套基因组上面竟然有两个定位！
查了好久的bug，终于搞清楚问题所在了！因为要对基因进行reads计数，所以要拿到基因在基因组上面的染色体起始终止坐标，结果发现了个十分诡异的现象，很多基因有多个坐标，比如下面这个PTPRS 在hg38这个基因组版本，居然有两个定位，因为我是写程序格式化得到的坐标，所以我check了我的程序，http://www.biotrainee.com/thread-472-1-1.html 感兴趣的同学可以点开看看我的代码！
十一
07
mysql的table居然有最大列限制
想着把TCGA的RPKM值矩阵表格写入到mysql，然后做一个查询网页给生物学家，我下载的是所有TCGA收集的mRNA表达数据集数据集-GSE62944 ，共9264个癌症样本，和741个正常组织的表达数据。当我想写入癌症表达矩阵的时候，报错了：
Error in .local(conn, statement, ...) :
could not run statement: Too many columns
简单搜索了一下，发现是mysql有最大列的限制，但是我不是很懂计算机，所以没太看明白该如何调整参数使得mysql列限制扩充：http://dev.mysql.com/doc/refman/5.7/en/column-count-limit.html 所以就把癌症表达矩阵根据癌症拆分了，癌种数量如下：
十一
07
【直播】我的基因组（五）:测试数据及参考基因组的准备
我的全基因组数据还没拿到，而且还会推迟，简单说(tu)明(cao)一下原因（还好当初为了避免广告嫌疑一直没说是哪个公司负责测序，反正用的是illumina的hiseqX10这个测序啦，所以可以尽情的吐槽）。
心烦意乱的吐槽线 Continue reading
十一
03
热图最佳实践-pheatmap
用pheatmap来绘图首先要安装这个包，它就一个功能，画出热图即可，号称是pretty heatmap，的确比其它的好用很多。我以前写过《一步一步学习heatmap.2》的教程，很简单的那种，所以就没有公布在博客上面，结果发现很多其它博客居然能先我一步发出。其实包括本次的pheatmap指南，都没什么好发，的在R里面也是傻瓜式出图，无法就是自己熟练一下参数而已，又不是开发一个包，没什么技术含量。我这里单独提一下pheatmap是因为它的确非常好用，将会是我画热图的不二之选。比如下面这个，是我最喜欢的： Continue reading