Daily Archives: 2015年5月21日

五 21

Biostrings包简介

Posted on by

首先讲讲它的对象

有下面几个字符串对象BString, DNAString, RNAString and AAString可以通过以下代码构造它们：

b <- BString("I am a BString object")

d <- DNAString("TTGAAAA-CTC-N")

这两个对象的区别是DNAstring对象对字符串的要求严格一些，只有IUPAC字符和+-字符可以。

对构造好的对象可以通过下标来取子字符串对象，也可以通过subseq来取，但是子字符串仍然是数据对象，只有通过toString函数才能把它们转化成字符串。

用length(dd2)和nchar(toString(dd2))都可以找到我们Biostrings对象的长度。但是后者速度会很慢。

Views(RNAString("AU"), start=0, end=2)这个函数可以把string对象任意截取成list

start, end and width可以作用于我们截取的list，判断list里面的元素在原来的string对象上面的起始终止及长度信息。

 

接下来讲这个包带有的一个比对函数！

> pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede")

Global PairwiseAlignmentsSingleSubject (1 of 2)pattern: [1] succ--eed subject: [1] supersede score: -33.99738

> pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede",type = "local")

Local PairwiseAlignmentsSingleSubject (1 of 2)pattern: [1] su subject: [1] su score: 5.578203

> pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede",gapOpening = 0, gapExtension = 1)

Global PairwiseAlignmentsSingleSubject (1 of 2)pattern: [1] su-cce--ed subject: [1] sup--

可以看出这个比对函数可以调整的参数实在是太多了，而且改变参数之后比对情况大不一样，还有很多参数就不一一细讲了。

这个比对结果可以赋值给一个变量，保存比对的对象。

psa1 <- pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede")

class(psa1)

summary(psa1)

class(pattern(psa1))

class(summary(psa1))

score(psa2)

还可以自己构建打分矩阵来进行比对。

submat <-

+ matrix(-1, nrow = 26, ncol = 26, dimnames = list(letters, letters))

diag(submat) <- 0

Biostrings包简介1454

psa2 <-pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede",substitutionMatrix = submat, gapOpening = 0, gapExtension = 1)

我们的包还自带了两个非常流行的氨基酸比对矩阵PAM和BLOSUM

ls("package:Biostrings")可以查看这个包所有的对象。

data(package="Biostrings")可以查看这个包所有的数据对象

还有很多其它函数

还可以去除adaptor，挺好玩的

既然有配对比对函数，那么就有多重比对函数！

我们可以读取clustaW, Phylip and Stolkholm这几种不同的比对结果文件来构造多重比对对象。

library(Biostrings)这个包里面自带了两个文件，我们可以示范一下构建对象。

origMAlign <- readDNAMultipleAlignment(filepath = system.file("extdata", "msx2_mRNA.aln", package="Biostrings"), format="clustal")

phylipMAlign <- readAAMultipleAlignment(filepath = system.file("extdata","Phylip.txt", package="Biostrings"),format="phylip")

 

对构造好的多重比对对象就可以构建进化树啦，代码如下！

sdist <- stringDist(as(origMAlign,"DNAStringSet"), method="hamming")

> clust <- hclust(sdist, method = "single")

> pdf(file="badTree.pdf")

> plot(clust)

> dev.off()

Biostrings包简介2345

五 21

Bioconductor的DO.db包介绍

Posted on by

Bioconductor的包都是同样的安装方法：

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

还有GO.bd包是完全一模一样的规则！！！

加载这个包可以发现它依赖于好几个其它的包，这也是我比较喜欢R的原因，它会自动把它需要的包全部安装加载进来，不需要自己一个个调试！

> library(DO.db)

载入需要的程辑包：AnnotationDbi

载入需要的程辑包：stats4

载入需要的程辑包：GenomeInfoDb

载入需要的程辑包：S4Vectors

载入需要的程辑包：IRanges

载入需要的程辑包：DBI

> help(DO.db)

> ls("package:DO.db")

[1] "DO"          "DO_dbconn"   "DO_dbfile"   "DO_dbInfo"   "DO_dbschema" "DOANCESTOR"  "DOCHILDREN"  "DOID"        "DOMAPCOUNTS"

[10] "DOOBSOLETE"  "DOOFFSPRING" "DOPARENTS"   "DOSYNONYM"   "DOTERM"      "DOTerms"     "Secondary"   "show"        "Synonym"

[19] "Term"

这个包里面有19个数据对象！都是比较高级的S4对象。

比如我们可以拿DOTERM[1:10]这个小的数据对象来做例子！example=DOTERM[1:10]

因为example是一个高级对象，所以无法直接查看，需要用as.list方法来查看

> as.list(example)

$`DOID:0001816`DOID: DOID:0001816Term: angiosarcomaSynonym: DOID:267Synonym: DOID:4508Synonym: "hemangiosarcoma" EXACT []Secondary: DOID:267Secondary: DOID:4508

~~~~~~~~~~~~共十个DO条目

对每一个DO条目来说都有DOID,Term,Synony这些函数可以取对应的值。

下面是对DO的有向无环图的数据解读

xx <- as.list(DOANCESTOR)可以查看每个DO与它所对应的上级条目DO，每个DO都会有不止一个的上级DO。

xx <- as.list(DOPARENTS)可以查看每个DO与它所对应的父条目DO，每个DO都有且只有一个父DO。

xx <- as.list(DOOFFSPRING)可以查看每个DO与它所对应的下级DO的关系列表，大多数DO都不止一个子条目DO，所有的下级DO都会列出。

xx <- as.list(DOCHILDREN)以查看每个DO与它所对应的子条目DO的关系列表，大多数DO都不止一个子条目DO。

还有Lkeys(DOTERM)可以查看数据库里面的所有的DO条目的ID号

> head(keys(DOTERM))

[1] "DOID:0000000" "DOID:0001816" "DOID:0002116" "DOID:0014667" "DOID:0050004" "DOID:0050012"

dbmeta(GO_dbconn(), "GOSOURCEDATE")

可以查看这个DO库的制备时间

> dbmeta(DO_dbconn(), "DOSOURCEDATE")

[1] "20140417"

五 21

RNA-seq流程对基因和转录本的表达量的计算

Posted on by

bedtools multicov和htseq-count都可以用来对基因和转录本的表达量的计算！！！

我们总共有四个样本，已经比对到小鼠的mm9基因组上面了，数据大小如下

RNA-seq流程对基因和转录本的表达量的计算111

然后对基因和转录本计数需要一些额外的信息，即各个基因及转录本的位置信息，gtf文件需要在UCSC等各大数据库下载

RNA-seq流程对基因和转录本的表达量的计算170

然后我们制作一个config文件配置我们的数据地址

cat sample_bam.config 可以看到文件内容如下

/data/mouse/ptan/740WT1.bam

/data/mouse/ptan/741WT2.bam

/data/mouse/ptan/742KO1.bam

/data/mouse/ptan/743KO2.bam

几个批处理文件名及内容分别如下

bedtools_multicov.sh  bedtools_multicov_transcript.sh  htseq.sh  htseq_transcript.sh

 

while read id

do

echo $id

new=`echo $id |cut -d"/" -f 5`

echo $new

bedtools multicov -bams $id -bed /data/mouse/mouse_mm9_gene.bed  > $new.gene.bedtools_multicov.count

done <$1

 

while read id

do

echo $id

new=`echo $id |cut -d"/" -f 5`

echo $new

bedtools multicov -bams $id -bed /data/mouse/mouse_mm9_transcript.bed  > $new.transcript.bedtools_multicov.count

done <$1

 

while read id

do

echo $id

new=`echo $id |cut -d"/" -f 5`

htseq-count -f bam $id /data/mouse/Mus_musculus.NCBIM37.67.gtf.chr  > $new.gene.htseq.count

done <$1

 

while read id

do

echo $id

new=`echo $id |cut -d"/" -f 5`

htseq-count -f bam --idattr transcript_id $id /data/mouse/Mus_musculus.NCBIM37.67.gtf.chr  > $new.transcript.htseq.count

done <$1

 

批量运行这些程序后就能对它们分别分情况进行计数，也能比较这两种计数方法的区别！

RNA-seq流程对基因和转录本的表达量的计算1201

 

可以看出区别还是很大的！！！

RNA-seq流程对基因和转录本的表达量的计算1219

我肯定没搞懂它们的原理，这完全就不一样，已经不是区别的问题了！！！

对于每个个体输出的计数文件，接下来就可以用DESeq等包来进行差异基因分析啦！

五 21

SAMStat软件使用说明书

Posted on by

这个软件是对我们的比对结果（通常是bwa,bowtie,tophat,hisat,star）bam或者sam来进行一个可视化的总结，类似于fastqc对我们的fastq测序结果做一个可视化总结，非常好用。

一．下载并安装该软件

软件主页是http://samstat.sourceforge.net/ 里面对该软件进行非常详细的说明

包括installation和usage，我这里简单的翻译一下。

Wget http://liquidtelecom.dl.sourceforge.net/project/samstat/samstat-1.5.tar.gz

解压开看里面的readme有介绍如何安装这个软件

Unpack the tarball:

bash-3.1$ tar -zxvf samstat-XXX.tar.gz

bash-3.1$ cd samstat

bash-3.1$ ./configure

bash-3.1$ make

bash-3.1$ make check

bash-3.1$ make install

如果用root命令就可以直接用samstat啦

如果没有root权限，安装的时候稍微有点不同

./configure  --prefix=/home/jmzeng/my-bin/

make

make install

很简单的

二，数据，就是我们的bam文件啦

三，运行命令

SAMStat软件使用说明书710

四，结果

简单看看samtools flagstat 740WT1.bam  的结果

19232378 + 0 in total (QC-passed reads + QC-failed reads)

0 + 0 duplicates

18846845 + 0 mapped (98.00%:-nan%)

0 + 0 paired in sequencing

0 + 0 read1

0 + 0 read2

0 + 0 properly paired (-nan%:-nan%)

0 + 0 with itself and mate mapped

0 + 0 singletons (-nan%:-nan%)

0 + 0 with mate mapped to a different chr

0 + 0 with mate mapped to a different chr (mapQ>=5)

然后再看看我们的samstat的结果！

740WT1.bam.samstat.html

一个网页，非常丰富的内容

SAMStat软件使用说明书1168

内容太多了，我懒得解释了

见软件说明书http://davetang.org/wiki/tiki-index.php?page=SAMStat

五 21

使用Bedtools对RNA-seq进行基因计数

Posted on by

以前是没有想过用这个软件的，直到有一个我的htseq无法对比对的bam文件进行基因计数（后来我才发现htseq无法计数的原因是gtf版本不同导致坐标不同，而且gtf对染色体编号没有加上chr），我简单搜索了一下，发现bedtools multicov也有类似的功能，所以我选择它来试试看！

首先注意它需要sort的bam文件及bam的index

bedtools multicov depends upon index BAM files in order to count the number of overlaps in each BAM file. As such, each BAM file should be position sorted (samtool sort aln.bam aln.sort) and indexed (samtools index aln.sort.bam) with either samtools or bamtools.

首先安装它：

wget https://github.com/arq5x/bedtools2/releases/download/v2.23.0/bedtools-2.23.0.tar.gz

解压开后

Make clean

Make all

然后就可以看到它的bin目录下全部是程序啦

Bedtools使用笔记639

命令很简单的

bedtools multicov [OPTIONS] -bams BAM1 BAM2 BAM3 ... BAMn -bed  <BED/GFF/VCF>

By default, multicov will report the count of alignments in each input BAM file that overlap.

例子：

bedtools multicov -bams aln1.bam aln2.bam aln3.bam -bed ivls-of-interest.bed

ivls-of-interest.bed这个文件是必须的，可能需要自己制作，其实用gtf文件也可以的

chr1 0   10000   ivl1

chr1 10000   20000   ivl2

chr1 20000   30000   ivl3

chr1 30000   40000   ivl4

输出结果前三列是坐标，第四列是基因名，跟我们的bed文件一样，只是最后三列是三个样本的计数，是添加上来的！

chr1 0       10000   ivl1    100 2234    0

chr1 10000   20000   ivl2    123 3245    1000

chr1 20000   30000   ivl3    213 2332    2034

chr1 30000   40000   ivl4    335 7654    0

 

同样是对gene的reads计数，bedtools的multicov工具与htseq-count的区别是

i'd guess it's due to the fact that htseq-count only reports one hit per aligned read assuming that read is aligned uniquely and does not overlap multiple features listed in your GTF. if an aligned read hits more than one feature in your GTF then it doesn't report that hit. bedtools gives you raw hits which includes every 1 hit for every intersection of every alignment with any features in the GTF no matter how many times it aligned or how many features it hit. you might think, "wow, htseq-count is dropping a lot of information". yes, it is! i've moved to using other tools to count hits to genes (RSEM/eXpress) since they disambiguate those ambiguous alignments and as a result you get counts for all of your aligned reads. in a genome with alternative splicing you lose too much data using htseq-count, in my opinion.

而且专门有个文献在讨论这个问题

http://barcwiki.wi.mit.edu/wiki/SOPs/rna-seq-diff-expressions

http://www.nature.com/nbt/journal/v31/n1/abs/nbt.2450.html

Differential analysis of gene regulation at transcript resolution with RNA-seq

下面我讲一个实际的例子

我的bam文件如下

Bedtools使用笔记2406

bedtools multicov -bams 740WT1.bam 741WT2.bam 742KO1.bam 743KO2.bam -bed mm9.bed

Bedtools使用笔记2491

得到的这个矩阵就可以去用DESeq包来进行差异分析啦！