09

对vcf突变数据与公开发表的进行比对

当我们对NGS数据call了snp之后一般会输出成vcf格式的数据,一行代表一个突变,例如
20      2451451 .       G       T       1939.77 .
AC=1;AF=0.500;AN=2;BaseQRankSum=-10.134;DP=239;Dels=0.00;FS=2.276;HaplotypeScore=0.0000;MLEAC=1;MLEAF=0.500;MQ=60.00;MQ0=0;MQRankSum=-0.258;QD=8.12;ReadPosRankSum=0.823;SOR=0.870
GT:AD:DP:GQ:PL  0/1:150,89:239:99:1968,0,3874
#前几列记录着该突变发生在第几号染色体以及该染色体的哪个坐标,我们的参考基因组在该位点是什么碱基,我们测到的突变成了什么碱基。
最后两列是测序深度以及正负测序深度,或者ref和allele的测序深度。
只有第8列是最复杂的,可以有高达几百个数据信息,取决于我们用什么样的软件来call的snp,以及call了snp之后用什么样的软件做的注释。
接下来我们还需要探究我们找到的突变是否在其它以及公开发表的数据库里面被找到过,所以可以下载非常多的公共数据库进行比对,我所见过的有一下一些,估计完全下载有0.5T
dbsnp144 (这个是ncbi提供的最权威的啦)
cgi69
ExAC.vcf.gz(这个是broadinstitute提供的外显子联盟)
Cosmic_v73.ann.vcf.gz (这个是癌症突变信息集)
finalTCGA.vcf.gz (TCGA计划也是癌症相关的)
icgc.vcf.gz
dbNSFP2.6vcf
SCLP.ann.vcf.gz
CCLE.ann.vcf.gz
ESP6500-SIv2.vcf.gz (Variants from the Exome Sequencing Project (ESP))
adni-sum
safs-sum.indel.vcf.gz
gonl.vcf.gz
ssm.vcf.gz
ssi.vcf.gz
uk10k.vcf.gz
1000g-ph3v5.gff.gz  (千人基因组计划)
gwasCatalog.gff.gz  \
phewascatalog.gff.gz  \
dbgap-gwas.gff.gz  \
interproDomain.gff.gz \
clinvar.gff.gz \
RegulomeDB.gff.gz \
CancerGAMAdb.gff.gz \
29

推荐5个生物信息学领域的教授

排名不分先后:

推荐宾夕法尼亚州立大学的一个教授Istvan Albert

他写了一本书是: https://www.biostarhandbook.com/
他还可以授予网上课程学位:http://www.personal.psu.edu/iua1/certificate.html
他还推荐了一本R语言书籍:http://onepager.togaware.com/

关注一下华盛顿大学医学院的教授Obi L. Griffith

他的主页:http://www.obigriffith.org/

他的一个比较出名的的贡献是 www.rnaseq.wiki
他在 Biostars bioinformatics forum 非常活跃
他的课程包括Molecular Basis of Cancer (BIO5288) and Genetics and Genomics of Disease (BIO5487) at Washington University School of Medicine.
I was a TA for Genome Analysis (MEDG505) and the bioinformatics section of Advanced Human Molecular Genetics (MEDG520) and a guest instructor for Cell Biology For Biomedical Engineering Graduate Students (APSC552), Cell and Organismal Biology (BIOL111) and Cell Biology (BIOL200) at UBC.

关注一下华盛顿大学医学院的教授Malachi Griffith

他的个人主页是:http://www.malachigriffith.org/index.htm

他的github主页是:https://github.com/malachig

WashU TGI Faculty page: Profile
Linked In: Profile
Twitter: Feed
Google Scholar: Citations
Research Gate: Profile
Scopus: Profile
Open Research ID: Profile
Github: Profile
BioStar: Profile
SeqAnswers: Profile
Code Academy: Profile
Iterative Genomics Consulting: Company website
Flickr: Photostream
www.dgidb.org
www.alexaplatform.org

关注一下麦吉尔大学的Pablo Cingolani教授

他是snpeff的作者

他的github是:https://github.com/pcingola
现就职于McGill University

推荐弗吉尼亚大学的stephen教授

他是个人主页:http://stephenturner.us/

他所有公开的ppt : https://speakerdeck.com/stephenturner
stephen教授我要重点提一下,因为他的教育资源特别多。
29

关于qsub和condor两种在集群上面提交任务的方式比对

毕竟不是计算机科班出身,很多计算机概念我其实并不清楚,很多时候对一个任务的解决只是凭着自己的悟性来模仿接触到的做法,大多数情况下并无可厚非,能完成任务即可,但是碰到多人协作的环境,往往就出了问题。

我一直都知道可以 cat /proc/cpuinfo 来查看cpu数量,用free -g 来查看内存的总量及消耗,用top来查看当前运行的任务情况。
以前都是三到五人的环境公用一个服务器,一般都有几十个cpu和几百G的内存,而且服务器还经常空闲着,所以不会面临资源的问题。
而最近接触的项目是多个小组公用服务器,所以面临计算资源的分配,接触了qsub和condor两种不同的任务提交模式。
qsub和condor其实都是用来运行一个脚本的,它与我们用bash或者sh来运行一个脚本的区别在于是否立即执行我们的脚本。
比如我写了一个脚本myscript.sh
我可以 bash myscript.sh 然后里面看到这个脚本被运行了
但是用 qsub myscript.sh, 就只是向我们的集群提交了一个任务而已,这个脚本什么时候运行,占用多少内存和cpu来运行,需要特定的设置好。
这里有个并行计算的ppt:http://www.marquette.edu/mugrid/bootcamp/2010Summer/Parallel.pdf 从计算机专业的角度来讲解了两者的区别。
其实只要看了这个https://wikis.nyu.edu/display/NYUHPC/Tutorial+-+Submitting+a+job+using+qsub 就能很容易明白qsub的用法,毕竟它只是一个命令而已,但是我们需要记住的是下面几个命令,因为经常用
PBS commands covered in this section
qsub
submit a job for execution
qstat
examine the status of a job (we have discussed what this status may be)
qhold
put a job on hold
qrls
release a job
qsig
send a signal to a job
qdel
delete a job

而condor其实很少见的。

An example Condor submission file is shown below:
Executable = ./myexecutable
Universe = vanilla
Log = condor.log
Output = out.log
Error = err.log
should_transfer_files = YES
when_to_transfer_output = ON_EXIT
transfer_input_files = mydata.txt
Save your script to a file with a name such as mysubmitfile.condor.
Once you have saved your submission script you can submit using the following command:
/dcs/condor/condor/bin/condor_submit mysubmitfile.condor

condor_q 可以用来查看任务提交情况

condor_rm 可以用来杀掉提交的任务。

任务提交模式还是挺有用的

25

RNA-seq流程需要进化啦!

Tophat 首次被发表已经是6年前

Cufflinks也是五年前的事情了

Star的比对速度是tophat的50倍,hisat更是star的1.2倍。

stringTie的组装速度是cufflinks的25倍,但是内存消耗却不到其一半。

Ballgown在差异分析方面比cuffdiff更高的特异性及准确性,且时间消耗不到cuffdiff的千分之一

Bowtie2+eXpress做质量控制优于tophat2+cufflinks和bowtie2+RSEM

Sailfish更是跳过了比对的步骤,直接进行kmer计数来做QC,特异性及准确性都还行,但是速度提高了25倍

kallisto同样不需要比对,速度比sailfish还要提高5倍!!!

参考:https://speakerdeck.com/stephenturner/rna-seq-qc-and-data-analysis-using-the-tuxedo-suite

25

Shell里面的各种括号的区别

[],[[]],(),(()),{},{{}},以及在前面加上$的区别,以及它们互相杂交组合的区别!!!

[[ ]] double brackets

(())Double parentheses

{{}}double curly brackets

我们必须要记住的是下面

[] 相当于test,作逻辑判断

$( ) 与` ` (反引号) 都是用来做命令替换用

${ } 吧... 它其实就是用来作变量替换用的啦

(())就是用来计算的,相当于expr函数。

参考:http://sayle.net/book/

http://tldp.org/LDP/abs/html/index.html

 

我们首先看看一对的括号

首先[]是用来逻辑判断的,必须有空格

if [ -f binom.py ]

then

echo 'binom.py exists'

fi

或者

nub=$((i%4))

#echo $nub

if [ $nub == 0 ];then

echo "we need to sleep 4 hours"

sleep 14000

fi

这个[]操作符等价于test函数

if test $1 -gt 0
then
echo "$1 number is positive"
fi

但是都必须有空格!!!

参考:http://www.freeos.com/guides/lsst/ch03sec02.html

关于shell的test操作符还有很多http://tldp.org/LDP/abs/html/fto.html

( ) 将command group 置于 sub-shell 去执行,也称 nested sub-shell。

{ } 则是在同一个 shell 内完成,也称为non-named command group。

补充一个: {} 还可以做变量扩展 {5..9}  或者 {abcd}e, 自己运行一下就知道效果啦

这两个差异很小,而且一般用不着,就不讲了。

那么这一对的括号加上了$符号后又变成了上面鬼东西呢?

当然,只有$( ) ${ }才是合法的。

在 bash shell 中,$( ) 与` ` (反引号) 都是用来做命令替换用(command substitution)的。

在操作上,用$( ) 或` ` 都无所谓,用$( )的优点是:

1, ` ` 很容易与' ' ( 单引号)搞混乱,尤其对初学者来说

2, 在多层次的复合替换中,` ` 须要额外的跳脱( \` )处理,而$( ) 则比较直观

再让我们看${ } 吧... 它其实就是用来作变量替换用的啦。

一般情况下,$var 与${var} 并没有啥不一样。

但是用${ } 会比较精确的界定变量名称的范围,比方说:

[code][/code]

$ A=B

$ echo $AB

还可以用来截取变量,这个就很多花样啦

# 是去掉左边(在鉴盘上# 在$ 之左边)

% 是去掉右边(在鉴盘上% 在$ 之右边)

单一符号是最小匹配﹔两个符号是最大匹配

 

然后我们看看两对的括号:

nub=$((i%4)) 等价于$nub=`expr $i % 1` ;

((i++)) 等价于$i=`expr $i + 1` ;

所以(())就是用来计算的,而且里面的变量不需要$来标记啦

(在 $(( )) 中的变量名称,可于其前面加$ 符号来替换,也可以不用)

在(())前面加上$只是为了把计算结果给保存而已。

而两个中括号和两个大括号都是不合法的!

 

24

用 GMAP/GSNAP软件进行RNA-seq的alignment

软件的解说ppt :http://www.mi.fu-berlin.de/wiki/pub/ABI/CompMethodsWS11/MHuska_GSNAP.pdf

软件的下载地址: http://research-pub.gene.com/gmap/
有研究者认为这个软件的比对效果要比tophat要好,虽然现在已经多出来了非常多的RNA-seq的alignment软件,我还是简单看看这个软件吧,它本来是2005就出来的一个专门比对低通量的est序列,叫GMAP,后来进化成了GSNAP
step1:下载安装GMAP/GSNAP
是一个标准的linux源码程序,安装之前一定要看readme  ,http://research-pub.gene.com/gmap/src/README
解压进去,然后源码安装三部曲,首先 ./configu  然后make 最后make install
会默认安装在 /usr/local/bin 下面,这里需要修改,因为你可能没有 /usr/local/bin 权限,安装到自己的目录,然后把它添加到环境变量!
step2 :准备数据
比对一般都只需要两个数据,一是索引好的参考基因组,另一个是需要比对的测序数据。
但是这个GSNAP,还需要对应的GTF注释文件。
首先需要参考基因组:虽然软件本身提供了一个hg19的参考基因组,并且已经索引好了Human genome, version hg19 (5.5 GB)(http://research-pub.gene.com/gmap/genomes/hg19.tar.gz) ,但是下载很慢,而且不是对所有版本的GSNAP都适用。所以我这里对我自己的参考基因组进行索引。
gmap_build -D ./ -d  my_hg19.fa
然后取ensemble下载hg19的gtf文件。
然后还需要把自己下载的gtf文件也构建索引,需要两个步骤
cat my_hg19.gtf |  ~/software/gmap-2011-10-16/util/gtf_splicesites > my_hg19.splicesites
cat  my_hg19.splicesites  |   iit_store -o my_hg19.gtf.index
然后拷贝需要比对的RNA-seq测序文件
step3: 运行程序
就是一步比对而已
gsnap
-D /home/jschnable/gsnap_indexes/
-d arabidopsisv10
--nthreads=50
-B 5
-s  /home/jschnable/gsnap_indexes/arabidopsisv10.iit
-n 2
-Q
--nofails
--format=sam temp.fastq
> results.sam
参数有点多,自己看看说明书吧http://qteller.com/RNAseq-analysis-recipe.pdf 讲的非常详细。
24

用freebayes来call snps

step1:,软件安装
wget http://clavius.bc.edu/~erik/freebayes/freebayes-5d5b8ac0.tar.gz
tar xzvf freebayes-5d5b8ac0.tar.gz
cd freebayes
make
一个小插曲,安装的过程报错:/bin/sh: 1: cmake: not found
所以我需要自己下载安装cmake,然后把cmake添加到环境变量

首先下载源码包http://www.cmake.org/cmake/resources/software.html

wget http://cmake.org/files/v3.3/cmake-3.3.2.tar.gz

 解压进去,然后源码安装三部曲,首先 ./configu  然后make 最后make install  

cmake 会默认安装在 /usr/local/bin 下面,这里需要修改,因为你可能没有 /usr/local/bin 权限,安装到自己的目录,然后把它添加到环境变量!

step2:准备数据

an alignment file (in BAM format)
a reference genome in (uncompressed) FASTA format.
正好我的服务器里面有很多
不过,该软件也可以出了一个测试数据集
wget http://bioinformatics.bc.edu/marthlab/download/gkno-cshl-2013/chr20.fa
wget ftp://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/data/NA12878/alignment/NA12878.chrom20.ILLUMINA.bwa.CEU.low_coverage.20121211.bam
wget ftp://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/data/NA12878/alignment/NA12878.chrom20.ILLUMINA.bwa.CEU.low_coverage.20121211.bam.bai

用这个代码就可以下载千人基因组计划的NA12878样本的第20号染色体相关数据啦

step3:运行命令
网站给出的实例是:
freebayes -f chr20.fa \
    NA12878.chrom20.ILLUMINA.bwa.CEU.low_coverage.20121211.bam >NA12878.chr20.freebayes.vcf

一般就用默认参数即可

 
step4:输出结果解读
没什么好解读的了,反正是vcf文件,都看烂了,就那些东西
不过该软件的作者倒是拿该软件与broad用GATK做出的NA12878样本的突变数据做了比较

 

24

broad_institute收集的癌症数据

肾上腺皮质 Adrenocortical carcinoma ACC 92 Browse Browse
膀胱,尿路上皮 Bladder urothelial carcinoma BLCA 412 Browse Browse
乳腺癌 Breast invasive carcinoma BRCA 1098 Browse Browse
子宫颈 Cervical and endocervical cancers CESC 307 Browse Browse
胆管癌 Cholangiocarcinoma CHOL 36 Browse Browse
结肠腺癌 Colon adenocarcinoma COAD 460 Browse Browse
大肠腺癌 Colorectal adenocarcinoma COADREAD 631 Browse Browse
淋巴肿瘤弥漫性大B细胞淋巴瘤 Lymphoid Neoplasm Diffuse Large B-cell Lymphoma DLBC 58 Browse Browse
食管 Esophageal carcinoma ESCA 185 Browse Browse
FFPE试点二期 FFPE Pilot Phase II FPPP 38 None Browse
胶质母细胞瘤 Glioblastoma multiforme GBM 613 Browse Browse
脑胶质瘤 Glioma GBMLGG 1129 Browse Browse
头颈部鳞状细胞癌 Head and Neck squamous cell carcinoma HNSC 528 Browse Browse
肾嫌色 Kidney Chromophobe KICH 113 Browse Browse
泛肾 Pan-kidney cohort (KICH+KIRC+KIRP) KIPAN 973 Browse Browse
肾透明细胞癌 Kidney renal clear cell carcinoma KIRC 537 Browse Browse
肾乳头细胞癌 Kidney renal papillary cell carcinoma KIRP 323 Browse Browse
急性髓系白血病 Acute Myeloid Leukemia LAML 200 Browse Browse
脑低级神经胶质瘤 Brain Lower Grade Glioma LGG 516 Browse Browse
肝癌 Liver hepatocellular carcinoma LIHC 377 Browse Browse
肺腺癌 Lung adenocarcinoma LUAD 585 Browse Browse
肺鳞状细胞癌 Lung squamous cell carcinoma LUSC 504 Browse Browse
间皮瘤 Mesothelioma MESO 87 Browse Browse
卵巢浆液性囊腺癌 Ovarian serous cystadenocarcinoma OV 602 Browse Browse
胰腺癌 Pancreatic adenocarcinoma PAAD 185 Browse Browse
嗜铬细胞瘤和副神经节瘤 Pheochromocytoma and Paraganglioma PCPG 179 Browse Browse
前列腺癌 Prostate adenocarcinoma PRAD 499 Browse Browse
直肠腺癌 Rectum adenocarcinoma READ 171 Browse Browse
肉瘤 Sarcoma SARC 260 Browse Browse
皮肤皮肤黑色素瘤 Skin Cutaneous Melanoma SKCM 470 Browse Browse
胃腺癌 Stomach adenocarcinoma STAD 443 Browse Browse
胃和食管癌 Stomach and Esophageal carcinoma STES 628 Browse Browse
睾丸生殖细胞肿瘤 Testicular Germ Cell Tumors TGCT 150 Browse Browse
甲状腺癌 Thyroid carcinoma THCA 503 Browse Browse
胸腺瘤 Thymoma THYM 124 Browse Browse
子宫内膜癌 Uterine Corpus Endometrial Carcinoma UCEC 560 Browse Browse
子宫癌肉瘤 Uterine Carcinosarcoma UCS 57 Browse Browse
葡萄膜黑色素瘤 Uveal Melanoma UVM 80 Browse Browse

看起来癌症很多呀,任重道远

24

perl的模块组织方式

如何使用自己写的私人模块

模块通俗来讲,就是一堆函数的集合。

Personally I prefer to keep my modules (those that I write for myself or for systems I can control) in a certain directory, and also to place them in a subdirectory. As in:

/www/modules/MyMods/Foo.pm
/www/modules/MyMods/Bar.pm

And then where I use them:

use lib qw(/www/modules);useMyMods::Foo; 

useMyMods::Bar;

As reported by "perldoc -f use":

It is exactly equivalent to
BEGIN { require Module; import Module LIST; }
except that Module must be a bareword.

Putting that another way, "use" is equivalent to:

  • running at compile time,
  • converting the package name to a file name,
  • require-ing that file name, and
  • import-ing that package.

So, instead of calling use, you can call require and import inside a BEGIN block:

BEGIN{require'../EPMS.pm';
  EPMS->import();}

And of course, if your module don't actually do any symbol exporting or other initialization when you call import, you can leave that line out:

BEGIN{require'../EPMS.pm';}
比如我的一个模块如下,命名为my_stat.pm
package my_stat;
sub mean{
my $sum=0;
$sum+=$_ foreach @_;
$sum/($#_+1);
}
#print &mean(1..10),"\n";
sub stddev{
$avg=&mean(@_);
#print "$avg\n";
my $sum=0;
$sum+=($_-$avg)**2 foreach @_;
sqrt($sum/($#_));
#It will be different if you use $#_+1;
#sqrt($sum/($#_+1));
}
#print &stddev(1..10),"\n";
1;
里面有我定义好的两个函数 mean 和 stddev , 那么我就可以在我的其它perl程序里面直接引用这个模块,从而使用我的两个自定义函数。
use lib "./";  #取决于你把自定义模块my_stat.pm放在哪个目录
use my_stat;
print my_stat::stddev(1..10),"\n";
18

一个基因坐标定位到具体基因的程序的改进

这是为了回答以前的一个疑问:任意给定基因组的 chr:pos, 判断它在哪个基因上面?这个程序难吗?

基因的chr,start,end都是已知的

学术一点讲述这个问题:已知CNV数据在染色体上的position如chr1:2075000-2930999,怎样批量获取其对应的Gene Symbol呢(批量)

数据如下:

head gene_position.hg19 //21629

1 chr19 58858171 58874214 A1BG ENSG00000121410

2 chr12 9220303 9268558 A2M ENSG00000175899

3 chr12 9381128 9386803 A2MP1 ENSG00000256069

9 chr8 18027970 18081198 NAT1 ENSG00000171428

10 chr8 18248754 18258723 NAT2 ENSG00000156006

12 chr14 95058394 95090390  ENSG00000273259

13 chr3 151531860 151546276 AADAC ENSG00000114771

14 chr2 219128851 219134893 AAMP ENSG00000127837

15 chr17 74449432 74466199 AANAT ENSG00000129673

16 chr16 70286296 70323412 AARS ENSG00000090861

head pfam.df.hg19.bed   //340960

chr1  12190          12689          Helicase_C_2     0        +        12190          12689          255,255,0

chr1  69157          69220          7tm_4        0        +        69157          69220          255,255,0

chr1  69184          69817          7TM_GPCR_Srsx         0        +        69184          69817          255,255,0

chr1  69190          69931          7tm_1        0        +        69190          69931          255,255,0

chr1  69490          69910          7tm_4        0        +        69490          69910          255,255,0

现在需要对我们的pfam数据进行注释,根据每一行的chrpos来看看是属于哪一个基因

总共会有338879 pfam记录可以注释上基因。

注释之后应该是 head pfam.gene.df.hg19  这个样子

CDK11B       chr1  1571423     1573930     Pkinase        0        -         1571423     1573930     255,255,0

CDK11B       chr1  1572048     1573921     Pkinase_Tyr          0        -         1572048     1573921     255,255,0

CDK11B       chr1  1572120     1572823     Kinase-like  0        -         1572120     1572823     255,255,0

CDK11B       chr1  1572120     1572820     Kinase-like  0        -         1572120     1572820     255,255,0

CDK11B       chr1  1572120     1572817     Kinase-like  0        -         1572120     1572817     255,255,0

CDK11B       chr1  1573173     1573918     Kinase-like  0        -         1573173     1573918     255,255,0

CDK11B       chr1  1575747     1577317     Daxx  0        -         1575747     1577317     255,255,0

CDK11B       chr1  1576417     1577347     Nop14         0        -         1576417     1577347     255,255,0

CDK11B       chr1  1576423     1577332     Mitofilin      0        -         1576423     1577332     255,255,0

CDK11B       chr1  1576432     1577317     SAPS  0        -         1576432     1577317     255,255,0

我的第一个程序用的是全基因位点扫描到hash的方法。这样需要扫描13,1390,4974个位点,多于三分之一的基因组,这样是非常浪费内存的,尤其是keys需要多个字节。

我用了256G的服务器都没有运行完。

后来我取巧了把我的gene_position.hg19文件用split命令分成了25个,然后循环25次对pfam.df.hg19.bed  文件进行注释。

 

这样的确可以解决问了,而且只需要32G的内存的服务器即可,时间也很快,就十多分钟吧。

但这只是取巧的方法,应该要从算法上面优化,首先我仅仅做一个改动,就是不再扫描全基因的位点,对每个基因,我以1K的窗口来取位点进行扫描。这样我判断pfam的坐标时候,也以1K为最小单位进行判断。

这样只需要不到30s就可以出结果,总共注释了303474pfam记录,还不是最终的338879,因为我这次只注释了基因的1000整数倍基因区间,这样如果pfam记录落在一个基因的起始终止点不到1K位置时就不会被注释。这时候需要对代码进行继续优化。

 脚步懒得上传了,在我的有道云笔记里面。

http://note.youdao.com/share/?id=58e66d138e9434284ffa61c53b65abdc&type=note

09

affymetix的基因表达芯片数据差异基因分析

我主要是看了一个差异分析的教程,讲的非常详细,全面,我先简单列出这个教程,然后再贴出我的代码

GEO本来只有三种层级的数据,分别是Sample, Platform, and Series
现在共有14,927 platforms,包括主流的affymetrix,agilent,illumina等产商的芯片,以及它们在不同领域的应用(snp,snv,gwas等等),以及各种不同的生物体(人,小鼠,大鼠)
这个分析流程,仅仅针对于affymetrix公司的基因表达相关的芯片数据。
目录如下:
因为他也是转载,所以链接失效了,现在的链接如下:
其实根据目录名重新搜索肯定能得到内容的, 链接失效太正常了。
具体内容,我整理并且重新注释了以下,在有道云笔记里面。
基本上只需要用心看这个教程,都能上手芯片数据的差异分析,但这只是差异分析的一种方法而已,而且还是非常过时的方法。
现在比较流行DESeq,edgeR等高通量测序的差异分析包,即使是十几年前的芯片数据,也不需要下载cel那种数据,可以直接下载每个项目的表达量矩阵Series Matrix File(s)
然后在R里面用read.table,调整好参数就可以直接读取啦!
07

liftover基因组版本直接的coordinate转换

下载地址:http://hgdownload.cse.ucsc.edu/admin/exe/
使用方法:【从hg38转到hg19】
因为主流的基因组版本还是hg19,但是时代在进步,已经有很多信息都是以hg38的形式公布出来的了。
比如,我下载了pfam.df这个protein domain注释文件,对人的hg38基因组每个坐标都做了domain注释,数据形式如下:
查看文件内容head pfam.hg38.df ,如下:
PFAMID chr start end strand
Helicase_C_2 chr1 12190 12689 +
7tm_4 chr1 69157 69220 +
7TM_GPCR_Srsx chr1 69184 69817 +
7tm_1 chr1 69190 69931 +
7tm_4 chr1 69490 69910 +
7tm_1 chr1 450816 451557 -
7tm_4 chr1 450837 451263 -
EPV_E5 chr1 450924 450936 -
7TM_GPCR_Srsx chr1 450927 451572 -
我想把domain的起始终止坐标转换成hg19的,就必须要借助UCSC的liftover这个工具啦
这个工具需要一个坐标注释文件 http://hgdownload-test.cse.ucsc.edu/goldenPath/hg38/liftOver/
而且它只能对bed等符合要求的格式进行转换
示例如下:
chr7  127471196  127472363  Pos1  0  +  127471196  127472363  255,0,0
chr7  127472363  127473530  Pos2  0  +  127472363  127473530  255,0,0
很简单的,把自己的文件随便凑几列信息,做成这个9列的格式即可
cat pfam.hg38.df |sed 's/\r//g' |awk '{print $2,$3,$4,$1,0,$5,$3,$4,"255,0,0"}'  >pfam.hg38.bed
这样就有了足够的文件可以进行坐标转换啦,转换的命令非常简单!
chmod 777 liftOver
 ./liftOver pfam.hg38.bed hg38ToHg19.over.chain pfam.hg19.bed unmap
然后运行成功了会有 提示,报错一般是你的格式不符合标准bed格式,自己删掉注释行等等不符合的信息即可
Reading liftover chains
Mapping coordinates
转换后,稍微检查一下就可以看到坐标的确发生了变化,当然,我们只需要看前面几列信息即可
grep -w p53 *bed
pfam.hg19.bed:chr11 44956439 44959858 p53-inducible11 0 - 44956439 44959858 255,0,0
pfam.hg19.bed:chr11 44956439 44959767 p53-inducible11 0 - 44956439 44959767 255,0,0
pfam.hg19.bed:chr2 669635 675557 p53-inducible11 0 - 669635 675557 255,0,0
pfam.hg19.bed:chr22 35660826 35660982 p53-inducible11 0 + 35660826 35660982 255,0,0
仔细看看坐标是不是变化啦!
pfam.hg38.bed:chr11 44934888 44938307 p53-inducible11 0 - 44934888 44938307 255,0,0
pfam.hg38.bed:chr11 44934888 44938216 p53-inducible11 0 - 44934888 44938216 255,0,0
pfam.hg38.bed:chr2 669635 675557 p53-inducible11 0 - 669635 675557 255,0,0
pfam.hg38.bed:chr22 35264833 35264989 p53-inducible11 0 + 35264833 35264989 255,0,0
其实R里面的bioconductor系列包也可以进行坐标转换 http://www.bioconductor.org/help/workflows/liftOver/
这个可以直接接着下载pfam.df数据库来做下去。更方便一点。
我的数据如下,需要自己创建成一个GRanges对象

1

library(GenomicRanges)
pfam.hg38 <- GRanges(seqnames=Rle(a[,2]),
               ranges=IRanges(a[,3], a[,4]),
               strand=a[,5])
2
这样就OK拉,虽然这只是一个很简陋的GRanges对象,但是这个GRanges对象可以通过R的liftover方法来转换坐标啦。
library(rtracklayer)
ch = import.chain("hg38ToHg19.over.chain")

pfam.hg19 = liftOver(pfam.hg38, ch)

pfam.hg19 =unlist(pfam.hg19)
再把这个转换好的pfam.hg19 写出即可

 

06

JQuery学习笔记

以后写这样的文章就直接用有道云笔记分享啦,这样可以节约这个免费的云服务器的空间。

jquery学习笔记第一弹:基础语法

http://note.youdao.com/share/?id=82021515144eb4820762e9fdbc686340&type=note

JQuery笔记第二弹:ppt效果操作

http://note.youdao.com/share/?id=08eb606b2084b9b0d8c9eb5ef72e3433&type=note

JQuery笔记第三弹:操作html元素

http://note.youdao.com/share/?id=fb8ff7deeb186adb82751838bf82cfbe&type=note

JQuery笔记第四弹:循环,遍历,判断等语句实现

http://note.youdao.com/share/?id=746ac6f1a801351f49d13cb3d7a335bf&type=note

JQuery笔记第五弹:Ajax实现

http://note.youdao.com/share/?id=0b2c6fb8c89e307ec79602e6d67e7c66&type=note

JQuery参考手册-函数大全

http://note.youdao.com/share/?id=2e926f98c9bd51b1192d309706f8c1ca&type=note

 

 

01

2012-LAD的三个亚型的不同生物学意义

文献名:Differential Pathogenesis of Lung Adenocarcinoma Subtypes Involving Sequence Mutations, Copy Number, Chromosomal Instability, and Methylation
Lung adenocarcinoma (LAD)的遗传变异度很大。
这个癌症可以分成三类:The LAD molecular subtypes (Bronchioid, Magnoid, and Squamoid)
然后我们在三个subtypes里面分析了以下四个特征,发现不同subtypes差异非常显著。
1、Gene mutation rates (EGFR, KRAS, STK11, TP53),
2、chromosomal instability,
3、regional copy number
4、genomewide DNA methylation
另外三个临床特征也是很显著。
1、Patient overall survival,
2、cisplatin plus vinorelbine therapy response
3、predicted gefitinib sensitivity
所以,我们的分类非常好,而且对临床非常有帮助。
对LAD的研究数据包括
1,DNA copy number
2,gene sequence mutation
3,DNA methylation
4,gene expression
即使是TP53这样的基因在LAD的突变率也才35%,所以我们的LAD应该更加细分,因为EGFR mutation and KRAS mutation这样的突变对治疗很有指导意义,细分更加有助于临床针对性治疗方案的选择。
我们选取了116个LAD样本的数据,分析了1,genome-wide gene expression,,2,genomewide DNA copy number, 3,genome-wide DNA methylation, 4,selected gene sequence mutations
得到的结论是:LAD molecular subtypes correlate with grossly distinct genomic alterations and patient therapy response
数据来源如下:
Gene expression --> Agilent 44 K microarrays.
DNA copy number --> Affymetrix 250 K Sty and SNP6 microarrays.
DNA methylation --> MSNP microarray assay.
DNA from EGFR, KRAS, STK11 and TP53 exons --> ABI sequencers
我们用的是R语言包 ConsensusClusterPlus根据gene expression 来对我们的LAD进行分类molecular subtypes
分类的基因有506个(the top 25% most variable genes, 3,045, using ConsensusClusterPlus),A nearest centroid subtype predictor utilizing 506 genes

这三类LAD的过表达基因参与不同的生物功能,
Bronchioid – excretion genes, asthma genes, and surfactants (SFTPB, SFTPC, SFTPD);
Magnoid – DNA repair genes, such as thymine-DNA glycosylase (TDG);
Squamoid – defense response genes, such as chemokine ligand 10 (CXCL10)
而且也对应不同的临床数据
Bronchioid had the most females, nonsmokers, early stage tumors, and low grade tumors, the greatest acinar content, the least necrosis, and the least invasion.
Squamoid had the most high grade tumors, the greatest solid content, and the lowest papillary content.
Magnoid had themost smokers and the heaviest smokers by pack years.
它们的基因突变pattern也有很大区别。
Bronchioid had the greatest EGFR mutation frequency
Magnoid had the greatest mutation frequencies in TP53, KRAS and STK11.
为了研究不同亚型癌症的突变模式的不同(genomewide mutation rates),我们同时又研究了a large set of rarely mutated genes (n = 623) from the Ding et al. cohort

结论:
Bronchioid subtype 更有可能受益于EGFR inhibitory therapy
Magnoid tumors also have severe genomic alterations including the greatest CIN, the most regional CN alterations, DNA hypermethylation, and the greatest genomewide mutation rate.
the Squamoid subtype displayed the fewest distinctive alterations that included only regional CN alterations

31

2013-science-3205tumors-12types-4-ways-find-291HCD

文献名:Comprehensive identification of  mutational cancer driver genes across 12  tumor types
本文比较了四种寻找癌症驱动基因的方法,并且得到了综合性的、可靠的291个HCDs 基因列表。
数据来源于3205个肿瘤样本,共涉及到12种癌症。
 Cancer Gene Census (CGC) 数据库里面已经有了接近500个cancer genes
 癌症基因组研究分析可以得到数以万计的somatic mutations,但是其中很少一部分才是驱动肿瘤发生,发展的突变。
 而且大多数driver genes的突变频率很低,又由于肿瘤的异质性,大量样本的研究是必须的。
 主流的四种找癌症驱动基因的方法如下:
 1、Most common methods identify genes that are mutated more frequently than expected from the background mutation rate (recurrence)
 2、Other methods - a bias towards the accumulation of functional mutations (FM bias)
 3、other methods exploit the tendency to sustain mutations in certain regions of the protein sequence (CLUST bias)
 4、other approaches exploit the overrepresentation of mutations in specific functional residues, such as phosphorylation sites (ACTIVE bias)
 它们的代表软件是MuSiC, OncodriveFM, OncodriveCLUST and ActiveDriver
 本文把这四种方法进行了比较,并且综合了它们的结果。
 In summary, we provide a very reliable list of 291 HCDs and a second one, of 144 CDs, more comprehensive but with an expectedly higher false-positives rate
 One hundred and sixty-five of these candidates are novel findings not included in the CGC.
 然后,作者对这291个HCDs基因进行了功能分析,其中,它们主要集中在以下五个生物功能
Chromatin remodeling,
mRNA processing,
Cell signaling/proliferation,
Cell adhesion,
DNA repair/Cell cycle
然后把四种方法综合得到的291个HCDs基因与Cancer Gene Census (CGC) 数据库里面已经有的接近500个cancer genes进行综合比较
 本文首次展示了综合多种癌症驱动基因寻找方法的可能性,这种综合是基于两个事实:
 1,各种方法找癌症驱动基因本来就没有金标准,所以综合多种方法,更comprehensive。
 2,综合多种方法能更好的比较评估所找到的癌症驱动基因的准确性。
31

2014-4742samples-21tumors-Cancer5000-set-254-genes

文献名: Discovery and saturation analysis of  cancer genes across 21 tumour types.
我们知道对一个癌症的多个样本进行研究,其实很少高达20%样本突变 most intermediate frequencies (2–20%),还有很多低频突变,因为研究样本不够,从而不被发现
我们从 4,742个tumor-normal pairs的外显子测序数据集研究了somatic point mutations,共21种癌症。
癌症基因可能集中于以下七个功能:
proliferation,
apoptosis,
genome stability,
chromatin regulation,
immune evasion,
RNA processing
protein homeostasis
我们用有放回的抽样方法对数据进行统计,得出结论:如果我们对某个癌症的研究样本高达500-6000个的话,可以发现更多的临床低频突变。
这篇文章是为了解决以下三个问题:
1、大规模的研究cancer就能达到鉴别出所有的cancer driver genes的程度吗?(Coverage of known cancer genes)
2,增大样本量是否会揭示很多cancer driver genes?(Analysis of novel candidate cancer genes)
3、我们距离对所有的cancer driver genes的完全认知还有多远?(Saturation analysis)
突变数据的分析流程是Broad’s stringent filtering and annotation pipeline
突变情况如下:
3,078,483 somatic single nucleotide variations(SSNVs),
77,270 small insertions and deletions (SINDELs)
29,837 somatic di-, tri- or oligonucleotide variations (DNVs, TNVs and ONVs, respectively)
an average of 672 per tumour–normal pair
包括:
540,831 missense,
207,144 synonymous,
46,264 nonsense,
33,637 splice-site
2,294,935 non-coding mutations
我们找驱动基因的方法是:
We used the most recent version of the MutSig suite of tools
which looks for three independent signals:
high mutational burden relative to background expectation,
accounting for heterogeneity;
clustering of mutations within the gene;
enrichment of mutations in evolutionarily conserved sites.
我们把以上MutSig的几个独立组件分析得到的p-value组合起来,判断驱动基因,我们即对每种癌症做了单独分析,同时也对这21种癌症做了综合分析。
我们找到的驱动基因的结果:
单独对各个癌症进行分析,可以总共找到334个基因,当然不同癌症找到的基因有交集。
These 334 pairs involve 224 distinct genes.
The number of genes detected per tumour type varied considerably (range of 1–58)
找到的驱动基因的个数差异主要取决于癌症种类的不同,然后,跟该癌症的样本量有关。
只有22种基因能在超过三种癌症里面都是被判定为驱动基因。
如果我们把21种癌症合并起来找驱动基因,可以找到114个,其中有30个是单独对各个癌症进行分析所找不到的,有80个在单独癌症分析可以找到。
所以单独对各个癌症进行分析找到的224个基因里面,有140个是合并癌症分析找不到的。其实画一个韦恩图就很好理解了。
对各个癌症进行分析,共21次分析,加上合并分析,共22次飞行,总共可以得到a Cancer5000 set containing 254 genes.
我们再严格分析一下254个基因在Cancer5000 set,得到219 distinct genes.叫做Cancer5000-S (for ‘stringent’) gene。
 Cancer Gene Census (CGC)组织的 (v65)版本包含着130个cancer genes driven by somatic point mutations,其中82个被我们这次统计分析发现啦。
 Four genes encode anti-proliferative proteins, in which loss-offunction mutations would be expected to contribute to oncogenesis.
 Sixadditionalgenesencode proteins thatare clearlyinvolved incell  proliferation: RHEB, RHOA, SOS1, ELF3, SGK1 and MYOCD.
 Five genes encode pro-apoptotic factors, in which loss-of-function mutations would be expected to promote oncogenesis
Six genes encode proteins related to genome stability.
Fivegenesareassociatedwithchromatinregulation
Three genes encode proteins whose loss is expected to help tumours evade immune attack
Three genes are associated with RNA processing and metabolism.
One gene, TRIM23, is involved in protein homeostasis.
Beyond these 33 genes, the set of 81 novel genes is likely to contain
additional true cancer genes.
有返回抽样方法是:An effective test is to perform ‘down-sampling’; that is, to study how the number of discoveries increases with sample size, by repeating the analysis on random subsets of samples of various smaller sizes.
饱和度分析结果: 还远未到饱和,不同突变频率的基因被发现的个数随着样本量的增大而增多的速度不同。
Genes mutated in 20% of tumours are approaching saturation;
those mutated at frequencies of 10–20% are still rising rapidly, but at a decreasing rate;
those at 5–10%  increasing linearly;
and those at ,5% are increasingly at an accelerating rate.
我们对样本量的要求是:突变背景高的癌症(如,黑色素瘤)需要的样品更多,而那些突变背景低的癌症(如成神经细胞瘤)需要近650个样本就可以很好的分析驱动基因了
Creating a reasonably comprehensive catalogue of candidate cancer genes mutated in 2% of patients will require between approximately 650 samples (for tumours with ,0.5 mutations per Mb, such
as neuroblastoma) to approximately 5,300 samples (for melanoma, with 12.9 mutations per Mb)
31

2015-MADGiC-identify-cancer-driver-gene

最新的一个寻找cancer 的driver gene的软件:
Cancer is thought to result from the accumulation of causal  somatic mutations throughout the lifetime of an individual.
这些cancer-driving mutations 主要影响三类基因: 1、oncogenes 2、tumor-suppressor genes 3, stablity geens
第一个突变是tumorigenesis ,随后的突变就 driver tumor progression
识别这些突变非常有利于了解gene function 和药物靶点设计
区分 driver genes 和 passenger  genes 能更好的利用各种数据库得到的海量突变信息
基于频率的区分方法 rely on an estimate of a background mutation rate which  represents the rate of random passenger mutations.
也就是文献(Ding et al., 2008).提出的方法,但它忽略了以下四点
1、mutation type (transition versus transversion)
2、nucleotide context(which base is at the mutation site
3、dinucleotide context (which bases are located at neighboring sites to the mutation),
4、expression level of the gene
然后有文献提出了以下三种改进
Sjoblom et al.(2006) account for nucleotide and dinucleotide context in searching  for drivers of breast and colorectal cancer.
MuSiC (Dees et al.,2012) accounts for mutation type and allows for sample-specific mutation rates;
Lawrence et al.(2013) (MutSigCV) also allow for the inclusion of gene-specific factors such as expression level and replication timing.
但是他们有个共同延续下来的的缺点,就是默认驱动基因的突变频率要高于背景突变频率。
实际上,除了突变频率,还有一些criteria也很重要, 所以有两个数据库SIFT (first reported by Ng and Henikoff (2001), later updated by Kumar et al. (2009)),  Polyphen (Adzhubei et al., 2010)  和MutationAssessor (Reva et al., 2011)
这两个数据库整合了 sequence context, position, and protein characteristics to assess a mutation’s  functional impact.
总结一下identity cancer driver genes的criteria
1、mutation frequency,
2、mutation type,
3、gene-specific features such as replication timing and expression level that are known to affect background rates of mutation,
4、mutation-specific scores that assess functional impact, and the spatial patterning of mutations that only becomes apparent when thousands of samples are considered.
以前的方法都只是部分涉及到上面的criteria
而我们提出了a unified empirical Bayesian Model-based Approach for identifying Driver Genes in Cancer (MADGiC) that utilizes each of these features.
31

2014-REVIEW-identifying driver mutation in sequenced cancer genome

somatic  mutations 含义很广,包括:SNVs,Indel,CNAs,SVs等
However, the resulting sequencing datasets are large and complex, obscuring the clinically important mutations in a background of errors, noise, and random mutations.
Cancer is driven largely by somatic mutations that accumulate in the genome over an individual’s lifetime, with additional contributions from epigenetic and transcriptomic alterations
低通量时代研究,成功例子: imatinib has been used to target cells expressing the BCR-ABL fusion gene  in chronic myeloid leukemia
gefitinib has been used to inhibit the epidermal growth factor receptor in lung cancer
但远远不够。
NGS的三大挑战:1,indentying somatic mutations,误差/ 肿瘤异质性 2,识别driver genes 3,确定由somatic mutations 改变的pathways和其它生物过程
误差来源:optical PCR duplicates, GC-bias, strand bias (where reads indicating a possible mutation only align to one strand of DNA) and alignment artifacts resulting from low complexity or repetitive regions in the genome.
most methods for somatic mutation detection address only a subset of the possible sources of error,call snp的软件众多
identifying driver mutations的三个要点:
1,identifying recurrent mutations;
2,predicting the functional impact of individual mutations;
3,assessing combinations of mutations using pathways, interaction networks, or statistical correlations.
三个要点分别衍生了大量的软件,它们的问题在于:
1,直接看突变频率的那些软件to determine whether the observed number of mutations in the gene is significantly greater than the number expected according to a background mutation rate (BMR).
BMR 实在是太难确定了,低了会导致很多假阳性,而高了,又错过很多真实的driver mutations,但是突变频率非常高的那些基因肯定是没有问题的,比如说TP53,无论什么样的算法都会认为它是driver gene
2,考虑突变对蛋白功能的影响评分的那些软件,引入了一些先验假设:
evolutionary conservation,
known protein domains,
non-random clustering of mutations,
protein structure,
3,pathways, interaction networks, and de novo approaches的那些软件:
pathway(KEGG,GO,GSEA) 4个limitations,首先,大多数 annotated gene sets 包含的基因数太多,而我们的突变基因占该gene set的比例远达不到统计显著性。
然后,pathway并不是独立的,各个pathway之间的联系更重要
接着,把基因分割成pathway这样的小单元,忽略了单元外的联系
最后,只关注已知的 pathways, or gene set
过去的五年见证了癌症基因组测序研究翻天覆地的变化,但是距离它真正的临床应用还有以下几个挑战:
首先,我们忽略了non-coding somatic mutations
其次,很多我们定义的癌症种类其实是a mixture of these subtypes
然后,哪些癌症是可以合并研究的
最后,不同的NGS数据如何综合研究,包括WGS,WES,RNA sequencing, DNA methylation, and chromatin modifications
对某些患者来说,癌症精准医学已经来临,但是对大部分病人来说,前面的路还很长。