10

RNA-seq比对软件HISAT说明书

取代bowtie+tophat进行RNA-seq比对

HISAT全称为Hierarchical Indexing for Spliced Alignment of Transcripts,由约翰霍普金斯大学开发。它取代Bowtie/TopHat程序,能够将RNA-Seq的读取与基因组进行快速比对。这项成果发表在3月9日的《Nature Methods》上。

HISAT利用大量FM索引,以覆盖整个基因组。以人类基因组为例,它需要48,000个索引,每个索引代表~64,000 bp的基因组区域。这些小的索引结合几种比对策略,实现了RNA-Seq读取的高效比对,特别是那些跨越多个外显子的读取。尽管它利用大量索引,但HISAT只需要4.3 GB的内存。这种应用程序支持任何规模的基因组,包括那些超过40亿个碱基的。

HISAT软件可从以下地址获取:http://ccb.jhu.edu/software/hisat/index.shtml。

首先,我们安装这个软件!

Wget http://ccb.jhu.edu/software/hisat/downloads/hisat-0.1.5-beta-source.zip

官网下载的是源码包,需要make一下,make之后目录下面就多了很多程序,绿色的那些都是,看起来是不是很眼熟呀!!!

哈哈,这完全就是bowtie的模拟版本!!!

HISAT取代bowtie+tophat进行RNA-seq比对1222

也可以从github里面下载,wget https://codeload.github.com/infphilo/hisat/zip/master

下载后直接解压即可使用啦。当然这个软件本身也有着详尽的说明书

http://ccb.jhu.edu/software/hisat/manual.shtml

然后就是准备数据,它跟tophat一样的功能。就是把用RNA-seq方法测序得到的fastq文件比对到参考基因组上面,所以就准这两个文件了哦

接下来是运行程序!

说明书上面写着分成两个步骤,构建索引和比对。

这个软件包模仿bowtie自带了一个example数据,而且它的说明书也是针对于那个example来的,我也简单运行一下。

$HISAT_HOME/hisat-build $HISAT_HOME/example/reference/22_20-21M.fa 22_20-21M_hisat

构建索引的命令如上,跟bowtie一样我修改了一下

/home/jmzeng/hoston/RNA-soft/hisat-0.1.5-beta/hisat-build 22_20-21M.fa  my_hisat_index

连日志都跟bowtie一模一样,哈哈,可以看到我们的这个参考fasta文件 22_20-21M.fa 就变成索引文件啦,索引还是很多的!

HISAT取代bowtie+tophat进行RNA-seq比对1871

然后就是比对咯,还是跟bowtie一样

$HISAT_HOME/hisat -x 22_20-21M_hisat -U $HISAT_HOME/example/reads/reads_1.fq -S eg1.sam

我的命令是

/home/jmzeng/hoston/RNA-soft/hisat-0.1.5-beta/hisat -x  my_hisat_index -U ../reads/reads_1.fq  -S reads1.sam

1000 reads; of these:

1000 (100.00%) were unpaired; of these:

0 (0.00%) aligned 0 times

1000 (100.00%) aligned exactly 1 time

0 (0.00%) aligned >1 times

100.00% overall alignment rate

哈哈,到这里。这个软件就运行完毕啦!!!是不是非常简单,只有你会用bowtie,这个就没有问题。当然啦,软件还是有很多细节是需要调整的。我下面就简单讲一个实际的例子哈!

首先,我用了1.5小时把4.6G的小鼠基因组构建了索引

/home/jmzeng/hoston/RNA-soft/hisat-0.1.5-beta/hisat-build  Mus_musculus.GRCm38.fa.fa mouse_hisat_index

HISAT取代bowtie+tophat进行RNA-seq比对2512

然后对我的四个测序文件进行比对。

for i in *fq

do

/home/jmzeng/hoston/RNA-soft/hisat-0.1.5-beta/hisat  -x  /home/jmzeng/hoston/mouse/mouse_hisat_index  \

-p 30 -U  $i.trimmed.single  -S ./hisat_out/${i%.*}.sam

done

它运行的速度的确要比tophat快好多,太可怕的速度!!!!至于是否多消耗了内存我就没有看了

4.6G的小鼠,5G的测序数据,我只用了五个核,居然十分钟就跑完了!

然后听群友说是因为没有加 --known-splicesite-infile <path>这个参数的原因,没有用gtf文件来指导我们的RNA数据的比对,这样是不对的!

需要用下面这个脚本把gtf文件处理一下,然后导入什么那个参数来指导RNA比对。

extract_splice_sites.py genes.gtf > splicesites.txt

但是我报错了,错误很奇怪,没解决,但是我换了个 extract_splice_sites.py  程序,就可以运行啦!之前是HISAT 0.1.5-beta release 2/25/2015里面的python程序,后来我换做了github里面的就可以啦!

/home/jmzeng/hoston/RNA-soft/hisat-master/extract_splice_sites.py Mus_musculus.GRCm38.79.gtf >mouse_splicesites.txt

HISAT取代bowtie+tophat进行RNA-seq比对3218

21192819 reads; of these:
21192819 (100.00%) were unpaired; of these:
14236834 (67.18%) aligned 0 times
5437800 (25.66%) aligned exactly 1 time
1518185 (7.16%) aligned >1 times

感觉没有变化,不知道为什么?

21192819 reads; of these:

21192819 (100.00%) were unpaired; of these:

14236838 (67.18%) aligned 0 times

5437793 (25.66%) aligned exactly 1 time

1518188 (7.16%) aligned >1 times

32.82% overall alignment rate

发表这个软件的文献本身也把这个软件跟其它软件做了详尽的对比

http://www.nature.com/nmeth/journal/v12/n4/full/nmeth.3317.html

Program Run time (min) Memory usage (GB)
Run times and memory usage for HISAT and other spliced aligners to align 109 million 101-bp RNA-seq reads from a lung fibroblast data set. We used three CPU cores to run the programs on a Mac Pro with a 3.7 GHz Quad-Core Intel Xeon E5 processor and 64 GB of RAM.
HISATx1 22.7 4.3
HISATx2 47.7 4.3
HISAT 26.7 4.3
STAR 25 28
STARx2 50.5 28
GSNAP 291.9 20.2
OLego 989.5 3.7
TopHat2 1,170 4.3

 

 

 

参考:http://www.plob.org/2015/03/20/8980.html

http://nextgenseek.com/2015/03/hisat-a-fast-and-memory-lean-rna-seq-aligner/

 

10

RNA-seq的比对软件star说明书

类似于tophat的软件

首先当然是下载软件啦!

两个地方可以下载,一个是谷歌code中心,被墙啦,另一个是github,我的最爱。

wget https://codeload.github.com/alexdobin/STAR/zip/master

2013-RNA-seq的star-类似于tophat的软件说明书217

解压即可使用啦,其中程序在bin目录下面,根据自己的平台调用即可!

然后doc里面还有个pdf的说明文档,写的非常清楚,我也是看着那个文档学的这个软件!

接下来就是准备数据啦!

既然是类似于tophat一样的比对软件,当然是准备参考基因组和测序数据咯,毫无悬念。

然后 该软件也给出了一些测试数据

ftp://ftp2.cshl.edu/gingeraslab/tracks/STARrelease/2.1.4/

然后就是运行程序的命令!

分为两个步骤:首先构建索引,然后比对即可,中间的参数根据具体需要可以细调!

构建索引时候,软件说明书给的例子是

The basic options to generate genome indices are as follows:
--runThreadN NumberOfThreads
--runMode genomeGenerate
--genomeDir /path/to/genomeDir
--genomeFastaFiles /path/to/genome/fasta1 /path/to/genome/fasta2 ...
--sjdbGTFfile /path/to/annotations.gtf
--sjdbOverhang ReadLength-1

我模仿了一下。对我从ensembl ftp里面下载的老鼠基因组构建了索引

/home/jmzeng/hoston/RNA-soft/STAR-master/bin/Linux_x86_64/STAR  \

--runThreadN 30 #我的服务器还比较大,可以使用30个CPU  \

--runMode genomeGenerate \

--genomeDir  /home/jmzeng/hoston/mouse/STAR-mouse #构建好的索引放在这个目录 \

--genomeFastaFiles  /home/jmzeng/hoston/mouse/Mus_musculus.GRCm38.fa.fa \

--sjdbGTFfile /home/jmzeng/hoston/mouse/Mus_musculus.GRCm38.79.gtf \

--sjdbOverhang 284   #我的测序数据长短不一,最长的是285bp

当然注释的地方你要删除掉才行呀,因为cpu用的比较多。

2013-RNA-seq的star-类似于tophat的软件说明书1302

算一算时间,对4.6G的小鼠基因组来说,半个小时算是非常快的了!Bowtie2的index要搞两个多小时。

然后就是比对咯。这也是很简单的,软件说明书给的例子是

The basic options to run a mapping job are as follows:
--runThreadN NumberOfThreads
--genomeDir /path/to/genomeDir
--readFilesIn /path/to/read1 [/path/to/read2]

我稍微理解了一下参数,然后写出了自己的命令。

fq=740WT1.fq.trimmed.single

mkdir  740WT1_star

/home/jmzeng/hoston/RNA-soft/STAR-master/bin/Linux_x86_64/STAR  \

--runThreadN 20  \

--genomeDir  /home/jmzeng/hoston/mouse/STAR-mouse   \

--readFilesIn  $fq \

--outFileNamePrefix  ./740WT1_star/740WT1

如果输出文件需要被cufflinks套装软件继续使用。就需要用一下参数

Cufflinks/Cuffdiff require spliced alignments with XS strand attribute, which STAR will generate with --outSAMstrandField intronMotif option.

还有--outSAMtype参数可以修改输出比对文件格式,可以是sam也可以是bam,可以是sort好的,也可以是不sort的。

最后是输出文件解读咯!

其实没什么好解读的,输出反正就是sam类似的比对文件咯,如果还有其它文件,需要自己好好解读说明书啦。我就不废话了!

 

值得一提的是,该程序提供了2次map的建议

The basic idea is to run 1st pass of STAR mapping with the usual parameters , then collect the junctions detected in the first pass, and use them as ”annotated” junctions for the 2nd pass mapping.

在对RNA-seq做snp-calling的时候可以用到,尤其是GATK官方还给出了教程,大家可以好好学习学习!

http://www.broadinstitute.org/gatk/guide/article?id=3891

 

 

08

GWAS研究现状及资源下载

 GWAS研究是非常火的,NHGIR还专门为它开辟了专栏来介绍,下面这个图片也是来自于NHGIR组织,是GWAS近年来发表文章的状况。

图片1

可以在该文章上面下载这个所有的数据

wget http://www.genome.gov/admin/gwascatalog.txt

截至目前为止。2015年5月8日21:08:34

这个文档有19603行的数据,但是只有2113篇pubmed文献,共涉及到七千多个基因

有293种杂志都发过GWAS的文章,总共有2113篇文献,发表关联分析突变位点最多的是这篇文献23251661 在PLoS One杂志上面,共 949个rs突变

杂志排序
cut -f 2,5 gwascatalog.txt |perl -alne '{$hash{$_}++}END{print "$_" foreach sort {$hash{$a} <=> $hash{$b}} keys %hash}' |cut -f 2 |perl -alne '{$hash{$_}++}END{print "$_\t$hash{$_}" foreach sort {$hash{$a} <=> $hash{$b}} keys %hash}'
Hum Genet 41
Am J Hum Genet 62
Mol Psychiatry 64
PLoS One 132
PLoS Genet 145
Hum Mol Genet 168
Nat Genet 397
文章的rs突变点排序
cut -f 2,5 gwascatalog.txt |perl -alne '{$hash{$_}++}END{print "$_ $hash{$_}" foreach sort {$hash{$a} <=> $hash{$b}} keys %hash}'
24324551 PLoS One 241
24097068 Nat Genet 245
24816252 Nat Genet 299
23382691 PLoS Genet 699
23251661 PLoS One 949

数据打开如下:

我取了表头和第一行数据,然后把它转置了,这样方便查看

Date Added to Catalog 10/22/2014
PUBMEDID 24528284
First Author Ji Y
Date 08/01/2014
Journal Br J Clin Pharmacol
Link http://www.ncbi.nlm.nih.gov/pubmed/24528284
Study Citalopram and escitalopram plasma drug and metabolite concentrations: genome-wide associations.
Disease/Trait Response to serotonin reuptake inhibitors in major depressive disorder (plasma drug and metabolite levels)
Initial Sample Description 300 European ancestry Escitalpram treated individuals, 130 European ancestry Citalopram treated individuals
Replication Sample Description NA
Region 17q25.3
Chr_id 17
Chr_pos 79831041
Reported Gene(s) CBX4
Mapped_gene CBX8 - CBX4
Upstream_gene_id 57332
Downstream_gene_id 8535
Snp_gene_ids
Upstream_gene_distance 33.93
Downstream_gene_distance 2.12
Strongest SNP-Risk Allele rs9747992-?
SNPs rs9747992
Merged 0
Snp_id_current 9747992
Context Intergenic
Intergenic 1
Risk Allele Frequency 0.086
p-Value 2.00E-07
Pvalue_mlog 6.698970004
p-Value (text) (S-DCT concentration)
OR or beta NR
95% CI (text) NR
Platform [SNPs passing QC] Illumina [7,537,437] (Imputed)
CNV N

 

 

上面这个文件是由tab键分割的,每一列的意义如下!

Note: The SNP data in the catalog has been mapped to dbSNP Build 142 and Genome Assembly,

GRCh38/hg37.p13.

DATE ADDED TO CATALOG: Date added to catalog

PUBMEDID: PubMed identification number

FIRST AUTHOR: Last name of first author

DATE: Publication date (online (epub) date if available)

JOURNAL: Abbreviated journal name

LINK: PubMed URL

STUDY: Title of paper (linked to PubMed abstract)

DISEASE/TRAIT: Disease or trait examined in study

INITIAL SAMPLE SIZE: Sample size for Stage 1 of GWAS

REPLICATION SAMPLE SIZE: Sample size for subsequent replication(s)

REGION: Cytogenetic region associated with rs number (NCBI)

CHR_ID: Chromosome number associated with rs number (NCBI)

CHR_POS: Chromosomal position associated with rs number (dbSNP Build 132,

NCBI)

REPORTED GENE (S): Gene(s) reported by author

MAPPED GENE(S): Gene(s) mapped to the strongest SNP (NCBI). If the SNP is

located within a gene, that gene is listed. If the SNP is intergenic, the upstream and

downstream genes are listed, separated by a hyphen. UPSTREAM_GENE_ID:

Entrez Gene ID for nearest upstream gene to rs number, if not within gene (NCBI)

DOWNSTREAM_GENE_ID: Entrez Gene ID for nearest downstream gene to rs

number, if not within gene (NCBI)

SNP_GENE_IDS: Entrez Gene ID, if rs number within gene; multiple genes

denotes overlapping transcripts (NCBI)

UPSTREAM_GENE_DISTANCE: distance in kb for nearest upstream gene to rs

number, if not within gene (NCBI)

DOWNSTREAM_GENE_DISTANCE: distance in kb for nearest downstream

gene to rs number, if not within gene (NCBI)

STRONGEST SNP-RISK ALLELE: SNP(s) most strongly associated with trait +

risk allele (? for unknown risk allele). May also refer to a haplotype.

SNPS: Strongest SNP; if a haplotype is reported above, may include more than one

rs number (multiple SNPs comprising the haplotype)

MERGED: denotes whether the SNP has been merged into a subsequent rs record

(0 = no; 1 = yes; NCBI)

SNP_ID_CURRENT: current rs number (will differ from strongest SNP when

merged = 1)

CONTEXT: SNP functional class (NCBI)

INTERGENIC: denotes whether SNP is in intergenic region (0 = no; 1 = yes;

NCBI)

RISK ALLELE FREQUENCY: Reported risk allele frequency associated with

strongest SNP

P-VALUE: Reported p-value for strongest SNP risk allele (linked to dbGaP

Association Browser)

PVALUE_MLOG: -log(p-value)

P-VALUE (TEXT): Information describing context of p-value (e.g. females,

smokers).

Note that p-values are rounded to 1 significant digit (for example, a published pvalue of 4.8 x 10-7 is rounded to 5 x 10-7).

OR or BETA: Reported odds ratio or beta-coefficient associated with strongest

SNP risk allele

95% CI (TEXT): Reported 95% confidence interval associated with strongest SNP

risk allele

PLATFORM (SNPS PASSING QC): Genotyping platform manufacturer used in

Stage 1; also includes notation of pooled DNA study design or imputation of

SNPs, where applicable

CNV: Study of copy number variation (yes/no)

Updated: January 13, 2015

 

 

 

 

05

Bioconductor的数据包library(biomaRt)简介

 

这是发布在bioconductor平台上面的一个数据库文件,可以通过R里面下载安装并使用,非常方便。其实在ensembl数据库里面也有一个biomart,我之前也讲过这个平台,非常好用,可以把任意的数据库之间的ID号进行转换。

为了更好的理解和掌握biomaRt,我们可以先通过在线资源来了解一下它的原型biomart (http://www.biomart.org)。 biomart是为生物科研提供数据服务的免费软件,它为数据下载提供打包方案。它有许多成功的应用实例,比如欧洲生物信息学中心(The European Bioinformatics Institute ,EBI)维护的Ensembl数据库(http://www.ensembl.org/)就使用biomart提供数据批量下载服务, 还有COSMIC, Uniprot, HGNC, Gramene, Wormbase以及dbSNP等。

这个就是一个R平台的biomart而已,但是非常好用!

> library(biomaRt)

> head(listMarts(), 3)

biomart                           version

1    ensembl      ENSEMBL GENES 79 (SANGER UK)

2        snp  ENSEMBL VARIATION 79 (SANGER UK)

3 regulation ENSEMBL REGULATION 79 (SANGER UK)

这是这个biomart最具有代表性的三个数据库,用listMarts()可以查看得知,它总共有58个数据库。

ensembl <-  useMart("ensembl", dataset = "hsapiens_gene_ensembl")

这是创建了人的ensembl数据库对象

> head(listFilters(ensembl), 3)

name     description

1 chromosome_name Chromosome name

2           start Gene Start (bp)

3             end   Gene End (bp)

可以看到对人的数据库ensembl来说,有多种字段可以来挑选自己感兴趣的东西,最常用的的当然是染色体号及起始终止坐标啦,用listFilters(ensembl),以查看得知,它总共有284中挑选感兴趣数据的方式。

既然 chromosome_name是其中一个挑选字段,那么我们就可以看看,是如何进行挑选的

用filterOptions(myFilter, ensembl)可以看到它挑选参数非常之多,远不止我们所认为的染色体号码。

染色体号一般是1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,X,Y,MT

还有一堆稀奇古怪的标志,LRG_101,LRG_102,LRG_103,LRG_104,因为我们组装好的人的标准基因组还有很多小的片段不被计入染色体中。

然后还可以看到人的ensembl数据库对象,有很多的属性,最常见的当然是基因ID和转录本ID和蛋白的ID号啦!

> head(listAttributes(ensembl), 3)

name           description

1       ensembl_gene_id       Ensembl Gene ID

2 ensembl_transcript_id Ensembl Transcript ID

3    ensembl_peptide_id    Ensembl Protein ID

用listAttributes(ensembl),,以查看得知,它总共有1166个ID号,太恐怖了,我实在是没有想到!

 

那么接下来我简单讲讲这个包的几个应用吧

首先是根据entrez ID号来找

ensembl <-  useMart("ensembl", dataset = "hsapiens_gene_ensembl")

这样就得到了人的信息,然后我探究以下两个基因的其它信息。

entrzID=c("672","1")

getBM(attributes=c("entrezgene","hgnc_symbol","ensembl_gene_id"), filters = "entrezgene", values =entrzID, mart=ensembl )

entrezgene hgnc_symbol ensembl_gene_id

1          1        A1BG ENSG00000121410

2        672       BRCA1 ENSG00000012048

3        672       BRCA1         LRG_292

 

其实这个函数很简单,就是根据自定义的entrzID这个变量来找到一些数据,数据的属性是我自己定义的entrezgene","hgnc_symbol","ensembl_gene_id",所以它就显示这个信息给我,在我之前弄好的人的数据库里面寻找!listAttributes(ensembl),,以查看得知,它总共有1166个ID号,就是说,你可以挑选你想要的基因的1166种信息,包罗万象!!!

其它功能也是很简单的啦,自己多看帮助文档!

 

从上面的操作来看,使用biomaRt只需要两步,1,指定mart数据库,2,使用getBM获得注释。但是首先,我们如何知道有哪些服务器,以及这些服务器上哪些数据库呢?其次,我们如何获阳getBM中attributes,filters的正确设置呢?

关于第一个问题,我们可以使用biomaRt中的listMarts以及listDatasets两个函数来解决。

> marts <- listMarts(); head(marts) #查看当前可用的数据源 ,总共有58个数据源。

> ensembl <- useMart("ensembl") #使用ensembl数据源

> datasets <- listDatasets(ensembl); datasets[1:10,] #查看ensembl中可用数据库,共有69个物种的数据库!

对于第二个问题,我们使用biomaRt中的listFilters以及listAttributes两个函数来解决。

> mart <- useMart("ensembl", "hsapiens_gene_ensembl")  #首先使用人的数据库

>listAttributes(ensembl) #,以查看得知,它总共有1166个ID号,就是说,人的数据库可供挑选的信息多达1166种。

> filters <- listFilters(mart); filters[grepl("entrez", filters[,1]),] #总共有284中挑选感兴趣数据的方式。

最后的问题是,biomaRt会被如何使用呢?我们做注释的时候,怎么就想到要使用biomaRt呢?因为在注释上,各种ID,symbol, name之间的转换都可以考虑使用biomaRt来做。更重要的是,biomaRt还会有很多SNP, alternative splicing, exon, intron, 5’utr, 3’utr等等信息。当然,只要能做也数据库并使用SQL访问的数据都可以使用biomaRt来获取。所以我们的思路可以更加发散一些。

 

 

 

05

Bioconductor的数据包library(org.Hs.eg.db)简介

 

这是发布在bioconductor平台上面的一个数据库文件,可以通过R里面下载安装并使用,非常方便。而且用的是数据库存储方式,所以搜索起来也是非常快速。

这个包里面有28个主流数据资料文件,这样我们可以用select函数根据我们自己的ID在这28个数据库里面随意转换自己想要的信息!!!

当然我本人是比较喜欢直接下载原文件,然后写脚本自己进行各种数据直接的转换。

首先我们加载这个数据包,可以看到这个数据包依赖于很多其它的包,如果是第一次安装。会耗时很长!

Bioconductor的数据包org.Hs.eg.db269

用这个函数,可以看到这个org.Hs.eg.db数据对象里面包含着各大主流数据库的数据,一般人都比较熟悉的entrez ID 和ensembl 数据库的ID。

keytypes(org.Hs.eg.db)

##  [1] "ENTREZID"     "PFAM"         "IPI"          "PROSITE"

##  [5] "ACCNUM"       "ALIAS"        "ENZYME"       "MAP"

##  [9] "PATH"         "PMID"         "REFSEQ"       "SYMBOL"

##  [13] "UNIGENE"      "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"

##  [17] "GENENAME"     "UNIPROT"      "GO"           "EVIDENCE"

##  [21] "ONTOLOGY"     "GOALL"        "EVIDENCEALL"  "ONTOLOGYALL"

##  [25] "OMIM"         "UCSCKG"

然后,我们用select函数,就可以把任意公共数据库的数据进行一一对应了。

ensids <- c("ENSG00000130720", "ENSG00000103257", "ENSG00000156414",

"ENSG00000144644", "ENSG00000159307", "ENSG00000144485")

cols <- c("SYMBOL", "GENENAME")

select(org.Hs.eg.db, keys=ensids, columns=cols, keytype="ENSEMBL")

比如说,我们有几个ensembl的基因ID号。然后我们想找它所对应的gene名和缩略词简称,就通过select函数来搞定即可!

Bioconductor的数据包org.Hs.eg.db1158

select(org.Hs.eg.db, keys="BRCA1", columns=c("ENSEMBL","UNIGENE","ENTREZID","CHR","GO","GENENAME"), keytype="SYMBOL")

这样得到了这个BRCA1基因的大部分信息,只是它的GO条目太多了,看得有点乱。

Bioconductor的数据包org.Hs.eg.db1318

 

 

 

05

Bioconductor简介

主页:http://www.bioconductor.org/

文字介绍我懒得写了,具体大家参考

http://www.bioconductor.org/about/

http://blog.csdn.net/shmilyringpull/article/details/8542607

这是一个R语言进行生信分析的流程发布平台,每个包都解决生信的一个流程问题。到目前为止2015年5月5日10:57:29已经有了1024个包,所以大家可以看到生信分析是一个海量的任务了。每一个包都有着详尽的说明文档及脚本代码,还附带着数据,非常容易弄懂,接下来我会花一个月的时间好好学习这些包!

这1024个虽然还是R语言的包,但是它的安装方式与常规的R语言包已经有所区别了。

需要用一下代码来安装

source("http://bioconductor.org/biocLite.R")biocLite()

biocLite(c("GenomicFeatures", "AnnotationDbi"))

也是非常easy的。

现在这个平台上面有1024个包,241个实验数据,917个数据库文件!!!

We are pleased to announce Bioconductor 3.1,

consisting of 1024 software packages,

241 experiment data packages,

and 917 up-to-date annotation packages.

在MOOC上面有很多关于这个的公开课

http://bioconductor.org/help/course-materials/

 

这里面有很多生信方向的分析流程,包括了我之前提到了snp-calling,RNA-seq,CHIP-seq等等,当然最主要的还是芯片数据的处理。

Workflows »

Common Bioconductor workflows include:

这些流程基本上涉及到了现在生物信息学的主流方向,所以基本上掌握了这些包,就是一个合格的生物信息学人才啦!

更重要的是它有着917个数据库文件,里面的信息分门别类,几乎可以算作是生物信息学的百科全书啦!

主要的数据库包括以下。

 

Package Description
AnnotationHub Ensembl, Encode, dbSNP, UCSC data objects
biomaRt Ensembl and other annotations
PSICQUIC Protein interactions
uniprot.ws Protein annotations
KEGGREST KEGG pathways
SRAdb Sequencing experiments.
rtracklayer genome tracks.
GEOquery Array and other data
ArrayExpress Array and other data

 

 

 

 

 

 

 

 

 

 

 

 

 

05

RNA-seq完整学习手册!

需耗时两个月!里面网盘资料如果过期了,请直接联系我1227278128,或者我的群201161227,所有的资源都可以在 http://pan.baidu.com/s/1jIvwRD8 此处找到

搜索可以得到非常多的流程,我这里简单分享一些,我以前搜索到的文献。

 

RNA-seq完整学习手册141

北大也有讲RNA-seq的原理

链接:http://pan.baidu.com/s/1kTmWmv9 密码:6yaz

甚至,我还有个华大的培训课程!!!这可是5天的培训教程哦,好像当初还花了五千多块钱的资料!!!

链接:http://pan.baidu.com/s/1nt5OV5B 密码:gyul

RNA-seq完整学习手册294

优酷也有视频,可以自己搜索看看

RNA-seq完整学习手册312

然后还有几个pipeline,就是生信的分析流程,即使你啥都不会,按照pipeline来也不是问题啦

export PATH=/share/software/bin:$PATH

bowtie2-build ./data/GRCh37_chr21.fa  chr21

tophat -p 1 -G ./data/genes.gtf -o P460.thout chr21 ./data/P460_R1.fq  ./data/P460_R2.fq

tophat -p 1 -G ./data/genes.gtf -o C460.thout chr21 ./data/C460_R1.fq  ./data/C460_R2.fq

cufflinks -p 1 -o P460.clout P460.thout/accepted_hits.bam

cufflinks -p 1 -o C460.clout C460.thout/accepted_hits.bam

samtools  view  -h  P460.thout/accepted_hits.bam  >  P460.thout/accepted_hits.sam

samtools  view  -h  C460.thout/accepted_hits.bam  >  C460.thout/accepted_hits.sam

echo ./P460.clout/transcripts.gtf > assemblies.txt

echo ./C460.clout/transcripts.gtf >> assemblies.txt

cuffmerge -p 1 -g ./data/genes.gtf -s ./data/GRCh37_chr21.fa  assemblies.txt

cuffdiff -p 1 -u merged_asm/merged.gtf  -b ./data/GRCh37_chr21.fa  -L P460,C460 -o P460-C460.diffout P460.thout/accepted_hits.bam C460.thout/accepted_hits.bam

samtools  index  P460.thout/accepted_hits.bam

samtools  index  C460.thout/accepted_hits.bam

 

和另外一个

#!/bin/bash

# Approx 75-80m to complete as a script

cd ~/RNA-seq

ls -l data

 

tophat --help

 

head -n 20 data/2cells_1.fastq

 

time tophat --solexa-quals \

-g 2 \

--library-type fr-unstranded \

-j annotation/Danio_rerio.Zv9.66.spliceSites\

-o tophat/ZV9_2cells \

genome/ZV9 \

data/2cells_1.fastq data/2cells_2.fastq                  # 17m30s

 

time tophat --solexa-quals \

-g 2 \

--library-type fr-unstranded \

-j annotation/Danio_rerio.Zv9.66.spliceSites\

-o tophat/ZV9_6h \

genome/ZV9 \

data/6h_1.fastq data/6h_2.fastq                          # 17m30s

 

samtools index tophat/ZV9_2cells/accepted_hits.bam

samtools index tophat/ZV9_6h/accepted_hits.bam

 

cufflinks --help

time cufflinks  -o cufflinks/ZV9_2cells_gff \

-G annotation/Danio_rerio.Zv9.66.gtf \

-b genome/Danio_rerio.Zv9.66.dna.fa \

-u \

--library-type fr-unstranded \

tophat/ZV9_2cells/accepted_hits.bam                  # 2m

 

 

time cufflinks  -o cufflinks/ZV9_6h_gff \

-G annotation/Danio_rerio.Zv9.66.gtf \

-b genome/Danio_rerio.Zv9.66.dna.fa \

-u \

--library-type fr-unstranded \

tophat/ZV9_6h/accepted_hits.bam                      # 2m

 

# guided assembly

time cufflinks  -o cufflinks/ZV9_2cells \

-g annotation/Danio_rerio.Zv9.66.gtf \

-b genome/Danio_rerio.Zv9.66.dna.fa \

-u \

--library-type fr-unstranded \

tophat/ZV9_2cells/accepted_hits.bam                  # 16m

 

 

time cufflinks  -o cufflinks/ZV9_6h \

-g annotation/Danio_rerio.Zv9.66.gtf \

-b genome/Danio_rerio.Zv9.66.dna.fa \

-u \

--library-type fr-unstranded \

tophat/ZV9_6h/accepted_hits.bam                      # 13m

 

 

time cuffdiff -o cuffdiff/ \

-L ZV9_2cells,ZV9_6h \

-T \

-b genome/Danio_rerio.Zv9.66.dna.fa \

-u \

--library-type fr-unstranded \

annotation/Danio_rerio.Zv9.66.gtf \

tophat/ZV9_2cells/accepted_hits.bam \

tophat/ZV9_6h/accepted_hits.bam                        # 7m

 

head -n 20 cuffdiff/gene_exp.diff

 

sort -t$'\t' -g -k 13 cuffdiff/gene_exp.diff \

> cuffdiff/gene_exp_qval.sorted.diff

 

head -n 20 cuffdiff/gene_exp_qval.sorted.diff

05

国外最出名的R语言大会-useR

这是2014年的会议报告以及ppt,但是好像很多ppt都是需要翻墙才能下载

http://user2014.stat.ucla.edu/#tutorials

Morning Tutorials Monday, 9:15

Room Presenter Title
Palisades Salon A+B Max Kuhn Applied Predictive Modeling in R
Palisades Salon C+F Winston Chang Interactive graphics with ggvis
Palisades Salon D+E Yihui Xie Dynamic Documents with R and knitr [Slides] [Examples]
Hermosa Romain Francois C++ and Rcpp11 for beginners [slides]
Venice Bob Muenchen Managing Data with R
Sproul-Landing building, 3rd floor Matt Dowle Introduction to data.table [Tutorial] [Talk]
Sproul-Landing building, 4th floor Virgilio Gomez Rubio Applied Spatial Data Analysis with R
Sproul-Landing building, 5th floor Martin Morgan Bioconductor

Afternoon Tutorials Monday, 14:00

Room Presenter Title
Palisades Salon A+B Hadley Wickham Data manipulation with dplyr
Palisades Salon C+F Garrett Grolemund Interactive data display with Shiny and R
Palisades Salon D+E Drew Schmidt Programming with Big Data in R
Hermosa S繪ren H繪jsgaard Graphical Models and Bayesian Networks with R
Venice John Nash Nonlinear parameter optimization and modeling in R [slides]
Sproul-Landing building, 3rd floor Dirk Eddelbuettel An Example-Driven Hands-on Introduction to Rcpp [slides]
Sproul-Landing building, 4th floor Ramnath Vaidyanathan Interactive Documents with R
Sproul-Landing building, 5th floor Thomas Petzoldt Simulating differential equation models in R

 

然后2015年的也要开始了,有兴趣的朋友可以June 30 - July 3, 2015
Aalborg, Denmark看看,有很多干货分享!

http://user2015.math.aau.dk/#BN

2015的内容如下

 

04

topGO简单使用

首先载入这个包

source("http://bioconductor.org/biocLite.R")

biocLite("topGO")

biocLite("ALL")

library(topGO)

library(ALL)

data(ALL)

data(geneList)

data(GOdatat)

这样就载入了很多变量, ls()查看如下

[1] "affyLib"      "ALL"          "geneList"     "topDiffGenes"

其中ALL这个数据我在另一篇日志里面重点介绍了一下。

然后我简单提一下"geneList"

head(geneList)

1095_s_at   1130_at   1196_at 1329_s_at 1340_s_at 1342_g_at

1.0000000 1.0000000 0.6223795 0.5412240 1.0000000 1.0000000

str(geneList) 是一个向量,有323个数字。

Named num [1:323] 1 1 0.622 0.541 1 ...

- attr(*, "names")= chr [1:323] "1095_s_at" "1130_at" "1196_at" "1329_s_at" ...

然后简单查询该包的安装地址和一些文件

system.file(package = "topGO")

[1] "C:/Program Files/R/R-3.1.1/library/topGO"

在这个目录下面可以找到文件"examples/geneid2go.map"

里面的内容格式如下,第一列是gene的ID号,一般是entrez ID ,第二列是该基因所对应的GO所有的条目,用逗号分隔。

068724 GO:0005488, GO:0003774, GO:0001539, GO:0006935, GO:0009288

119608 GO:0005634, GO:0030528, GO:0006355, GO:0045449, GO:0003677, GO:0007275

此处省略一万行。

readMappings(file = system.file("examples/geneid2go.map", package = "topGO"))

这个函数可以读取我们的文件,返回一个list。是gene到go的映射,每个基因都有一个或者多个go条目。

这个list可以用inverseList这个函数反转,变成每个go条目到基因的映射。

构建topGO这个大全,需要的数据包括:

  1. 基因identifier,可以附上某种分数以便后面施用某种统计处理,分数可以是t检验的p值或者与某个表型的correlation等;
  2. identifier和GO term间的map,如果是芯片数据的话BioC里有多种注释包,声明包的名称即可。至于我等蛋白界苦人,也能自己构建map,见下;
  3. GO的层级结构,由GO.db提供,目前这个包只支持GO.db提供的结构:GOslim就再说了。

topGOdata对象构建函数的参数包括:

  1. ontology,可指定要分析的GO term的类型,即BP、CC之类;
  2. description:topGOdata的描述,可选;
  3. allGenes:基因identifier的原始列表,和后面的geneSelectionFun联合作用,得出参与分析的基因,可以是numeric或factor。
  4. geneSelectionFun:基因选择函数,如果前面allGenes是numeric的话就必须得指明此参数;
  5. nodeSize:被认为富集的GO term辖下基因的最小数目(>=),默认为1。
  6. annotationFun:基因identifier map到GO term的函数。

代码如下

BPterms <- ls(GOBPTerm)

geneID2GO=readMappings(file = system.file("examples/geneid2go.map", package = "topGO"))

直接使用系统自带的data(GOdata)数据,自己构建太麻烦了!

其实就是这就对ALL这个数据集来进行分析!!!

GOdata包含实例topGOdata对象。它可以用来直接运行富集分析。

topGOdata对象构建好后,即可利用这个包里的各种方法和函数做分析。

numGenes(GOdata) 查看对象包含的基因的数目

[1] 318

> description(GOdata)

[1] "Simple topGOdata object"

genes(GOdata)可以得到该对象里面所有的318个基因

geneScore() 可以得到想318个基因的分数

函数sigGenes()返回一个character vector,为各显著变化基因identifier。函数numSigGenes()则用于查看显著变化基因的数目。

函数updateGenes()可以修改topGOdata对象里所包含的基因。

想要看全部基因都对应上了哪些GO term,可用函数usedGO()得到一个character

 

基因集富集分析(gene set enrichment analysis)

首先看看GSEA的三种方式:

  1. 基于count,即仅要求输入一组基因,此种方式最为流行,可用Fisher's exact test、Hypegeometric  test和binomial test进行检验;
  2. 基于基因的score或rank,可用Kolmogorov-Smirnov like tests(即05年那篇PNAS的GSEA文章所用方法),Gentleman's Category、t-test等方法;
  3. 基于基因的表达,可从expression matrix直接分析,如Goeman's globaltest,以及GlobalAncova。

topGO提供两种分析方法,一种自由度更高但上手不易,本菜鸟还是跟着第二种走吧,较为用户友好但集成度较高。简单来说,就是用runTest()这个函数,要求三个主要的argument,一个是之前构建好的topGOdata类实例,第二个参数algorithm用于指定生成GO graph的方法,而参数statistic用于指定所用的检验方法,比如:

> resultFis <- runTest(GOdata, algorithm = "classic", statistic = "fisher")

> resultWeight <- runTest(GOdata, algorithm = "weight", statistic = "fisher")

> resultKS <- runTest(GOdata, algorithm = "classic", statistic = "ks")

> resultKS.elim <- runTest(GOdata, algorithm = "elim", statistic = "ks.elim")

runTest这一锤子买卖敲定后就能开始解读和展示结果了,得到的结果是topGOresult类的一个实例,其组成很简单,为对象的基本信息,以及各基因的分数(p值或者其他统计参数

 

 

我这里随便挑一个富集结果来看看

resultFis <- runTest(GOdata, algorithm = "classic", statistic = "fisher")

 

-- Classic Algorithm --

 

the algorithm is scoring 590 nontrivial nodes

parameters:

test statistic:  fisher

 

resultWeight <- runTest(GOdata, algorithm = "weight", statistic = "fisher")

 

-- Weight Algorithm --

 

The algorithm is scoring 590 nontrivial nodes

parameters:

test statistic:  fisher : ratio

然后我们对这两种富集方式来画图

pvalFis=score(resultFis) 得到矫正的P值

pvalWeight <- score(resultWeight , whichGO = names(pvalFis))

返回resultFis和resultWeight共有的基因在resultWeight中的分数。有了这两组分数,可以做一些比较,比如关联分析:

cor(pvalFis, pvalWeight)

[1] 0.370151

library(lattice)

xyplot(pvalWeight ~ pvalFis) 画了一个散点图

 

04

R语言里面的一个数据集ALL(Acute Lymphoblastic Leukemia)简介

这个数据内容太多了,我感觉自己也理解的不是很清楚!

非常多的R的bioconductor包都是拿这个数据集来举例子的,所以我简单的介绍一下这个数据集。

这个数据集是对ALL这个病的研究数据,共涉及到了128个ALL病人,其中95个是B细胞的ALL,剩余33个是T细胞的ALL。

是一个芯片数据,同时还包含有其它的病人信息。

大家首先要在R里面安装这个数据集

source("http://bioconductor.org/biocLite.R")

biocLite("ALL")

library(ALL)

data(ALL)

data(geneList)

在R里面输入str(ALL)可以看到这个数据的具体信息,但是非常多!

ALL

ExpressionSet (storageMode: lockedEnvironment)

assayData: 12625 features, 128 samples 

element names: exprs

protocolData: none

phenoData

sampleNames: 01005 01010 ... LAL4 (128 total)

varLabels: cod diagnosis ... date last seen (21 total)

varMetadata: labelDescription

featureData: none

experimentData: use 'experimentData(object)'

 pubMedIds: 14684422 16243790 

Annotation: hgu95av2

我们首先它的BT变量记录的是什么

R语言里面的一个数据集ALL750

可以看出它记录的是数据病人的分组信息。

bcell = grep("^B", as.character(ALL$BT))通过这句话可以挑选出B细胞病人

然后我们看看它的ALL$mol.biol变量记录是是什么

R语言里面的一个数据集ALL857

可以看到它记录的是这些病人的几种突变情况(molecular biology testing)

types = c("NEG", "BCR/ABL")

moltyp = which(as.character(ALL$mol.biol) %in% types)

用这个命令就能挑选出我们想研究的两组突变的病人。

然后我们还可以把刚才的两个标准用来从ALL数据集里面取想要的子集

ALL_bcrneg = ALL[, intersect(bcell, moltyp)]

 

 

同理我们可以查看这个数据集的非常多的变量信息。

包括sex,age,cod,diagnosis,等等,这个'data.frame':共有128 obs. of  21 variables:

R语言里面的一个数据集ALL1190

 

我们除了可以查看这个ALL数据集自带的变量,还可以通过一些方法来访问它的其它信息。

Exprs这个方法可以查看它的表达数据,可以看到有128个变量,12625行的探针数据。

str(exprs(ALL))

num [1:12625, 1:128] 7.6 5.05 3.9 5.9 5.93 ...

- attr(*, "dimnames")=List of 2

..$ : chr [1:12625] "1000_at" "1001_at" "1002_f_at" "1003_s_at" ...

..$ : chr [1:128] "01005" "01010" "03002" "04006" ...

 

还有很多很多函数都可以操作这个数据集,这样可以得到非常多的信息!我就不一一列举了

对这个数据的一系列操作可以画热图,见下面的教程!!!

http://www2.warwick.ac.uk/fac/sci/moac/people/students/peter_cock/r/heatmap/

 

03

生物信息入门视频推荐-新一代测序数据分析

  1. 新一代测序数据分析-在优酷里面可以搜索到,一下是配套视频的讲义及下载地址!
  2. Lecture Notes
  3. Lectures will appear below as they are presented. Homeworks are specified in each handout.
  4. Lecture 1 - slideshandouts. course information, homework and project information, introduction to computing, setting up you computer, basic unix command line usage, organizing your projects, homework 1.
  5. Lecture 2 - slideshandouts, The GFF formatsequence ontologies, basic Unix commands: wc, grep, cut, sort, redirecting input and output streams, piping commands, processing a tabular file with UNIX tools, homework 2
  6. Lecture 3 - slideshandouts. programming languages, download and install an proper editor, introduction to the AWK programming language, tabular file processing, filtering by feature types, Awk onliners explained, another collections of AWK oneliners, homework 3.
  7. Lecture 4 - slideshandouts, sequencing technologies, sequence representations, the FASTA format, processing FASTA files at the command line, homework 4.
  8. Lecture 5 - slideshandouts, string matching, edit distances, regular expressions, local and global alignments, homework 5.
  9. Lecture 6 - slideshandouts, introduction to using blast, legacy blast and blast+, preparing blast databases, performing a blastn query, formatting blast output, homework 6.
  10. Lecture 7 - slideshandouts, using blast, formatting databases, using the blastdbcmd, extract sequences, batch operations, formatting blast queries, homework 7.
  11. Lecture 8 - slideshandouts, blast score and E-values, search strategies, usage examples for blastn, blastp, blastx, tblastn, and tblastx, homework 8.
  12. Lecture 9 - slideshandouts, quality encodings, phred scales, the FASTQ format, homework 9.
  13. Lecture 10 - slideshandouts, file compression, gzip, zip, bz2, file archives, tarbombs, plotting fastq qualities homework 10.
  14. Lecture 11 - slideshandouts installing tools, quality control, adapter trimming, error corrections
  15. Lecture 12 - slideshandouts paired end sequencing, quality control for paired end sequencing, the bioawk language
  16. Lecture 13 - slideshandouts paired end sequencing, read stiching, automating tasks with shell scripts
  17. Lecture 14 - slideshandouts short read alignments, bwa, bowtie and other tools.
  18. Lecture 15 - slideshandouts the sequence alignment map SAM format
  19. Lecture 16 - slideshandouts the SAM/BAM format, sorting and indexing BAM files, using the samtools program
  20. Lecture 17 - slideshandouts aligning paired end reads, comparing and evaluating aligners, simulating sequencing reads with the wgsim tool
  21. Lecture 18 - slideshandouts read duplication, visualizing alignments with IGV and IGB
  22. Lecture 19, guest lecture by Nicholas Stoler - slides, the variant call format (VCF), calling variants with samtools mpileup
  23. Lecture 20,- slideshandouts origins of genome variations, more on SNP calling, successes and failures
  24. Lecture 21,- slideshandouts interval representation, BED and GFF formats, representing data
  25. Lecture 22,- slideshandouts interval operations: complement, extension, flanking, Using the BedTools package
  26. Lecture 23,- slideshandouts interval operations: intersect, window, selecting closest features
  27. Lecture 24,- slideshandouts an introduction to genome assembly, using the velvet assembler, evaluating genome assemblies with QUAST
  28. Lecture 25,- slideshandoutsmeta.tar.gz (25MB) an introduction to metagenomics, software packages mothur, QIIME and MetaSim, online tools RDP, MG-RAST
  29. Lecture 26,- slideshandoutslec26.tar.gz (25MB) an introduction to Chip-Seq technology, peak calling concepts, preprocessing and peak calling methods (part 1)
  30. Lecture 27,- slideshandouts, Chip-Seq peak calling sofware, preprocessing and peak calling methods (part 2)
  31. Lecture 28,- slideshandoutslec28.tar.gz basic RNA-Seq data analysis concepts, split read mapping
  32. Lecture 29, slideshandoutslec29.tar.gz RNA-Seq (part 2)
  33. Lecture 30, slideshandouts, bioinformatics beyond the command line: using R for data analysis
  34. Final Project 30, final-project, data for final project pony.tar.gz (17Mb) BMB 597D: Final project, 50% of the final grade, due 5pm Saturday Dec 14th, 2013
01

脚本作业-解读NCBI的ftp里面关于人的一些基因信息

为了感谢大家对我博客的关注,我在这里发布一个作业,适合菜鸟做的。里面有十几个类似的问题,大家可以下载数据自行处理,如果是问这些问题,我优先回答!

NCBI的ftp里面关于人的一些基因信息

我在NCBI的ftp服务器里面下载了这些数据,时间是2015年,大多是hg19系列的,文件名如下:

CDS.fa 这个是ensembl中人的CDS碱基序列文件,hg38

entrez2go.gene 这个是有go注释的基因情况,有一万八的基因都有go注释

entrez2name.gene 这个是NCBI的entrez ID号对应着基因名的文件

entrez2pubmed.gene 这个是NCBI的entrez ID号对应着该基因发表过的文章的ID号

entrez2refseq2ensembl.gene 这个是NCBI的entrez ID号对应着基因名的refseq的ID号和ensembl数据库的ID号

human_gene_info这个是基因的详细信息,包括基因的起始终止点坐标等等

Protein.fa 这个是ensembl中人的蛋白的氨基酸序列文件,有十万多个蛋白hg38

ref2ensembl.txt  这个是基因名的refseq的ID号和ensembl数据库的ID号

自行去NCBI的ftp服务器里面下载这些数据。

然后好好熟悉这些数据信息,回答一下几个问题:

人总的基因有多少个,它们分别分布在哪些染色体上面,基因的转录本分布情况如何,基因的长度分布如何,基因的外显子个数如何。

CD分子的基因有多少个,它们分别分布在哪些染色体上面,基因的转录本分布情况如何,基因的长度分布如何,基因的外显子个数如何。它们有没有氨基酸偏好性??

MHC系列基因信息?CCL系列基因信息如何?CXCL系列信息如何?或者你感兴趣的基因家族信息?

现在研究最热门的基因是什么?发表文章最多的前十个基因是什么?

基因长度情况如何?最长的基因多长?最短的基因多少bp,可靠吗?

蛋白质长度情况如何?

每条染色体的基因分别情况?基因在染色体那个地方分别最多?

请用图形展示你的结论!!!

 

如果你能回答以上问题,证明你的脚本水平不错了。

如果找不到我,看旁边的公告,加入生信菜鸟群,我就在里面!!!

01

CHIP-seq第三讲之使用MACS软件寻找peaks

在使用Bowtie比对于完Chip-Seq的结果后,就需要用到MACS或者ERANGE来找出峰所在的位置了。但是由于ERANGE的设置比较复杂,所以最为流行的还是MACS。

一.首先安装MACS软件

MACS有两个版本,分别是MACS14和MACS2。MACS2在很多方面都对MACS14做了重大改进,但目前还在测试阶段。我们依然以MACS14为例进行说明。

MACS软件的下载地址在wget https://codeload.github.com/taoliu/MACS/zip/master

这是一个python软件,有152M,已经算是很大了!所以需要按照安装python的方法来安装它!但是,好像这个是最新版的,我们还是用1.4版本吧

wget http://github.com/downloads/taoliu/MACS/MACS-1.4.2-1.tar.gz

其实它的readme已经把这个软件的各种安装使用方法讲的很清楚了。

https://github.com/taoliu/MACS/blob/master/README.rst

MACS软件的具体原理,大家去看文献,或者参考这篇文章

http://www.plob.org/2014/05/08/7227.html

很简单的一个python命令即可安装该软件python setup.py install --user

CHIP-seq第三讲之使用MACS软件寻找peaks752

二.然后准备该软件所需要的数据

是我们在前两篇文章中提到的数据

CHIP-seq第三讲之使用MACS软件寻找peaks786

三.接着运行MACS的命令

/home/jmzeng/.local/bin/macs14 -t Xu_WT_rep2_BAF155.fastq.trimmed.single.bam \

> -c Xu_WT_rep2_Input.fastq.trimmed.single.bam \

> -f BAM -g hs --bw 300 -w -S -n Xu_WT_rep2

CHIP-seq第三讲之使用MACS软件寻找peaks974

四.最后解读一下结果

CHIP-seq第三讲之使用MACS软件寻找peaks987

56K Apr 30 21:54 Xu_WT_rep2_model.r

5.5K Apr 30 22:21 Xu_WT_rep2_negative_peaks.xls

783K Apr 30 22:21 Xu_WT_rep2_peaks.bed

865K Apr 30 22:21 Xu_WT_rep2_peaks.xls

766K Apr 30 22:21 Xu_WT_rep2_summits.bed

唉,反正这也不是我的课题,懒得解释这些结果啦,等后来有机会再慢慢玩吧

 

 

参考 http://www.plob.org/2014/05/08/7227.html

附录:我们现在来了解如何设置参数。

参考自 http://www.plob.org/2014/01/26/7118.html

-t TFILE, –treatment=TFILE 输入文件名

-c CFILE, –control=CFILE 输入阴对文件名

-n NAME, –name=NAME 输入出文件名前缀

-f FORMAT, –format=FORMAT 输入文件格式,默认值为AUTO,可选的值为”BEG”,”ELAND”,”ELANDMULTI”,”ELANDMULTIPER”,”ELANDEXPORT”,”SAM”,”BAM”,”BOWTIE”等。

-g GSIZE, –gsize=GSIZE 比对模板大小。格式可以是:1.0e+9,或者1000000000,也可以缩写:’hs’ for 人类 (2.7e9), ‘mm’ for 大鼠(1.87e9), ‘ce’ for 线虫 (9e7) and ‘dm’ for 果蝇 (1.2e8), 默认值:hs

-s TSIZE, –tsize=TSIZE 设置为短序列的长度,默认值为25

-p PVALE, –pvalue=PVALUE 非峰可能性截取值,默认值为1e-5,这个值不能大太,超过0.9的话,可能无法输出正确的结果

-m MFOLD, –mfold=MFOLD 峰值高度相对于本底的比值,默认值为10,30。也就是说,最低值不能少于10,但比值超过30也不认为它是正常的一个峰。一般而言,低值设置为10是一个很好的区分点。如果这个值还是无法得到满意的结果,那么可以设置得更低,但最好还是使用–nomodel参数,使–nomodel设置为True,然后再传递–shiftsize及–bw参数给MACS。–shiftsize默认值为100,而–bw的默认值为300。

–diag 生成完整报表,会包括是否为真峰的可能性,但会严重拖累运算速度。

 

30

人的CD分子基因信息简介!

CD分子吧,它是Clusters of Differentiation的简写,是指一组分化抗原的家族,目前该家族已经有CD1——CD350甚至更多的成员.他们分布于T细胞等免疫细胞表面,参与免疫细胞各种表达,其中有整合素、受体、配体等蛋白分子,在免疫应答反应中参与识别、粘附和信号转导等功能.

我这里简单讲讲如何整理它们的基因信息,首先从NCBI里面下载的人的gene_info文件,然后通过脚本来查找CD分子信息。

perl -alne '{if (/\tCD\d+/ or /CD\d+\|/ ) {print}}' human_gene_info >CD.info

cut -f 2-5 CD.info >CD.table

再根据CD分子的排序把我们的信息重新排序

perl -alne '{/CD(\d+\w)/;$hash{$1}=$_}END{print $hash{$_} foreach sort {$a <=> $b}keys %hash}' CD.table >CD.table.sort

然后我发现了一个很有趣的问题,它们都是负义链上面的基因!

 

entrez ID gene symbol 正负链
911 CD1C - BDCA1|CD1|CD1A|R7
913 CD1E - CD1A|R2
909 CD1A - CD1|FCB6|HTA1|R4|T6
912 CD1D - CD1A|R3
910 CD1B - CD1|CD1A|R1
9266 CYTH2 - ARNO|CTS18|CTS18.1|PSCD2|PSCD2L|SEC7L|Sec7p-L|Sec7p-like
30011 SH3KBP1 - CD2BP3|CIN85|GIG10|HSB-1|HSB1|MIG18
23607 CD2AP - CMS
89886 SLAMF9 - CD2F-10|CD2F10|CD84-H1|CD84H1|SF2001
10849 CD3EAP - ASE-1|ASE1|CAST|PAF49
445347 TARP - CD3G|TCRG|TCRGC1|TCRGC2
915 CD3D - CD3-DELTA|IMD19|T3D
920 CD4 - CD4mut
922 CD5L - AIM|API6|PRO229|SP-ALPHA|Spalpha
925 CD8A - CD8|Leu2|MAL|p32
927 CD8BP - CD8B2
54675 CRLS1 - C20orf155|CLS|CLS1|GCD10|dJ967N21.6
3681 ITGAD - ADB2|CD11D
3683 ITGAL - CD11A|LFA-1|LFA1A
3684 ITGAM - CD11B|CR3A|MAC-1|MAC1A|MO1A|SLEB6
3687 ITGAX - CD11C|SLEB6
290 ANPEP - APN|CD13|GP150|LAP1|P150|PEPN
115708 TRMT61A - C14orf172|GCD14|Gcd14p|TRM61|hTRM61
2526 FUT4 - CD15|ELFT|FCT3A|FUC-TIV|FUTIV|LeX|SSEA-1
2215 FCGR3B - CD16|CD16b|FCG3|FCGR3|FCR-10|FCRIII|FCRIIIb
4055 LTBR - CD18|D12S370|LT-BETA-R|TNF-R-III|TNFCR|TNFR-RP|TNFR2-RP|TNFR3|TNFRSF3
930 CD19 - B4|CVID3

 

 

 

 

 

 

 

30

自学CHIP-seq第二讲之过滤数据并比对

这个是有着非常成熟的流程了,我就不细讲了!

我们随机挑选两个文件来跑一下CHIP-seq的流程吧,其中一个是.部分进行免疫共沉淀前的DNA(input DNA)作为空白对照。

5.5G Apr 30 10:31 Xu_WT_rep2_BAF155.fastq

18G Feb 13 20:37 Xu_WT_rep2_Input.fastq

首先进行质量控制,过滤低质量的reads

这里我选取的是DynamicTrim.pl 和

脚本如下

for id in *fastq

do

echo $id

perl DynamicTrim.pl $id

done

接下来

for id in *.trimmed

do

echo $id

perl LengthSort.pl $id

Done

这样就得到了过滤后的reads,可以进行比对啦!

图片1

当然,中间文件可以删掉啦,不然太占空间了,我还只是取了两个数据,要是把这个文章的八个数据都跑完就太纠结了。

然后用bowtie比对

#samtools faidx hg19.fa

#Bowtie2-build hg19.fa hg19

for i in *single

do

bowtie2 -x /home/jmzeng/ref-database/hg19 -U $i -S  $i.sam

samtools view -bS $i.sam> $i.bam

done

输出的bam文件就需要用MASC这个软件来找peak了

30

自学CHIP-seq第一讲之文献解读

我这里选择的CHIP-seq文章题目是

CARM1 Methylates Chromatin Remodeling Factor BAF155 to Enhance Tumor Progression and Metastasis

文章链接http://www.sciencedirect.com/science/article/pii/S1535610813005369

这是2013年的文章,算是蛮新的了,主要探究了CARM1这个基因

然后我简单搜索了一些这个基因的信息

9606 10498 CARM1

- PRMT4

MIM:603934|HGNC:HGNC:23393|

Ensembl:ENSG00000142453|HPRD:09158|Vega:OTTHUMG00000180699

19 19p13.2 coactivator-associated arginine methyltransferase 1

protein-coding CARM1 coactivator-associated arginine methyltransferase histone-arginine methyltransferase CARM1|protein arginine N-methyltransferase 4 20150308

该基因是多种肿瘤相关的转录因子的共激活剂(激活蛋白;转录辅助激活蛋白;转录共同活化子)。

文章作者做了以下四件事

Knockout of CARM1 Using ZFN in Breast Cancer Cells

Identification of BAF155 as a Novel CARM1 Substrate

Methylation of BAF155 Promotes Tumor Growth and Metastasis

Methylated BAF155 Gains Unique Chromatin Association

 

所以就有两种细胞,一种是野生型WT,一种是突变的MUT细胞

然后它们分别做了两个重复,一种是input一种是BAF155免疫测序。

CHIP-seq一定是有一个input对照文件,和一个真正的免疫共沉淀的测序文件。

这样就有八个测序文件。

我们随机挑选两个文件来跑一下CHIP-seq的流程吧,其中一个是.部分进行免疫共沉淀前的DNA(input DNA)作为空白对照。

5.5G Apr 30 10:31 Xu_WT_rep2_BAF155.fastq

18G Feb 13 20:37 Xu_WT_rep2_Input.fastq

然后我随便在网上找了一个生信分析流程

  1. 标准信息分析
    a)   对测序数据进行base calling、raw data 数据整理及数据质量评估;
    b)   去接头污染,去低质量reads和产量情况统计
    c)   Bisulfite 测序序列与参考基因组序列的比对
    d)   深度和覆盖度分析
    e)   C 碱基的甲基化水平
    f)   全基因组甲基化水平分布趋势
  2. g)  全基因组DNA甲基化图谱
  3. h)  差异性甲基化区域(DMR)分析

 

参考

http://www.plob.org/2012/09/29/3760.html

http://www.plob.org/2012/01/09/1605.html

http://www.plob.org/2012/01/08/1538.html

 

30

阿里巴巴免费的服务器体验好差!

不知道为什么最近进入自己的网页后台总是很慢,发个日志也慢,很是郁闷!

本来以为是免费的空间快用完了,所以慢,结果一查,根本就没有用多,其实我很想投诉一下阿里巴巴!

想想该搞个国外服务器了,然后把网站搬家!

QQ截图20150430150101

30

Figtree的把进化树文件可视化

下载软件

http://tree.bio.ed.ac.uk/software/figtree/

我们这里就在window平台下使用,所以直接下载zip包即可,解压即可使用

准备数据

我这里就简单的用muscle生成了一个树文件来看看结果TRAV.fa 是一百多个类似的基因

muscle -in TRAV.fa -out tmp

muscle -maketree -in tmp  -out TRAV.tree

这个树文件TRAV.tree是Newick format,可以直接被figtree识别从而画图

软件使用

很简单,下载,点击即可使用,然后导入树文件,就可以直接出图啦!

Figtree的把进化树文件可视化368

30

Muscle进行多序列比对

软件的主页是

http://www.drive5.com/muscle/

进入主页,简单看看软件介绍,这个软件还是蛮牛的,一个人在家里自己写出来的,当然,对于普通人来说,这个软件跟clustalW没什么区别,反正都是多序列比对啦!

我们下载适合我们平台的版本即可!

Muscle进行多序列比对193

准备数据,我这里选择的是几个短小的蛋白

Muscle进行多序列比对215

 

这里有两种比对方式,都是很简单的命令

一种是先比对,再生成树文件(树的格式是Newick format, )

muscle -in mouse_J.pro -out mouse_J.pro.a

muscle -maketree -in mouse_J.pro.a -out mouse_J.phy (这里有两种构建树的方式)

另外一种是比对成aln格式的数据,然后用其它软件(phyml或者phylip)来生出树文件

muscle -in mouse_J.pro   -clwout seqs.aln

可以看到比对的效果还是蛮好的,是aln格式的比对文件,这个格式非常常用

Muscle进行多序列比对505

或者输出phy格式的比对文件,

muscle -in mouse_J.pro  -physout seqs.phy

Muscle进行多序列比对685

可以被phyml等软件识别,然后来构建进化树,见  http://www.bio-info-trainee.com/?p=626

21

美国Minnesota大学的生信全套课件分享

刚才在知乎什么看到了一篇分享pacbio的数据特征,顺便看到了Minnesota大学的关于生物信息的教程的ppt合集,所以就想打包下载。

https://www.msi.umn.edu/tutorial-materials

这个网页里面有64篇pdf格式的ppt,还有几个压缩包,本来是准备写爬虫来爬去的,但是后来想了想有点麻烦,而且还不一定会看,反正也是玩玩
就用linux的命令行简单实现了这个爬虫功能。
curl https://www.msi.umn.edu/tutorial-materials >tmp.txt
perl -alne '{/(https.*?pdf)/;print $1 if $1}' tmp.txt >pdf.address
perl -alne '{/(https.*?txt)/;print $1 if $1}' tmp.txt
perl -alne '{/(https.*?zip)/;print $1 if $1}' tmp.txt >zip.address
wget -i pdf.address
wget -i pdf.zip
这样就可以啦!
教程ppt列表如下,大家有兴趣的可以自行下载浏览。

2009-04-22-mrm-presentation_0.pdf               Matlab_viz_image_UMR.pdf
Analyzing ChIP at the command line.pdf          MaxQuant_Introduction_112409.pdf
Analyzing ChIP using Galaxy.pdf                 Maxquant-step-by-step_rs091124.pdf
Badalamenti_PacBio_tutorial_12-10-2014.pdf      MSI Applications Catalog Oct 21 MB slides.pdf
basics_chip_seq.pdf                             MSIIntro2013Jun18.pdf
Best_Practices_GATK_Variant_Detection_v1_0.pdf  MSIIntroBMEN5311.pdf
blast2go.pdf                                    MSI_Workshop_for_Introduction_to_Structure_based_Drug_Design.pdf
ClinProTools_0.pdf                              MTLB_GPUs.pdf
CUDA_Programming.pdf                            OpenMP.tutorial_1.pdf
cuda_tutorial_performance.pdf                   Open_Source_Proteomics_1.pdf
FLUENT_2009April21_final.pdf                    OptimizingWithGA.pdf
FLUENT_tutorial_2008aug14fin.pdf                Orbi_Data_Analysis_092811.pdf
galaxy_101_V4_ljm_0.pdf                         Partek Training Handout_miRNA and mRNA Data Analysis.pdf
GPU_tools.pdf                                   PerformanceTuning_itasca_11_27_12_0.pdf
gpututorial-msi.pdf                             PETSc_Tutorial.pdf
Hands_On_Tutorial_Using_ProTIP.pdf              Phi_Intro.pdf
Introduction to MSI Systems.pdf                 Protein_Grouping_FDR_Analysis_and_Database_Pratik_March2012_Draft.pdf
Introduction_to_PEAKS_0.pdf                     Proteomics_MSI_072309_Print.pdf
Introduction_to_SBDD.pdf                        pymol_v5.pdf
IntroMPI2011july19c.pdf                         QC_illumina_galaxy_V1_ljm.pdf
IntroMPI2012_July25-part1.pdf                   Quality Control of Illumina Data at the Command Line.pdf
IntroMSI2014.pdf                                remotevisualization.pdf
IntroNWChem.pdf                                 RISS_Hsapiens_variant_Detection_v3.0-small.pdf
IntroOpenMP_2011jun28b.pdf                      RNA_seq_Lecture2_2014_v2.pdf
Intro_to_GAMESS.pdf                             RNA-Seq mod1v6.pdf
IntroToGaussian09.pdf                           R_Spring2012_ver2.pdf
introtomolpro.pdf                               SchrodingerTutorial2011.pdf
Intro_to_MSI_Physicists.pdf                     Sybyl.pdf
intro-to-perl.pdf                               Tutorial-Hsap-v15.pdf
Matlab_11_29_UMR.pdf                            Tutorial-Stuber-v12-1.pdf
Matlab_PCT.pdf                                  unix2013.6.18.pdf
MATLAB_Tuning.pdf                               WRKSP_2_19.pdf

Total wall clock time: 40m 22s
Downloaded: 64 files, 249M in 40m 2s (106 KB/s)

我都已经下载好了,打包压缩到群里面啦!