生物信息学生R入门教程–输入与输出

Bioconductor应用,读写测序文件

现在有一个sequence文件(比如fasta, fastq, BAM, gff, bed, 或者 wig)需要读入R,应该怎么办?对于序列文件,我们需要Bioconductor的Biostrings包。对于BAM等文件,Rsamtools以及rtracklayer软件包。

> library("Biostrings")
> ## 读写fasta文件
> filepath <- system.file("extdata", "someORF.fa", package = "Biostrings")
> fasta.info(filepath, seqtype = "DNA")
## YAL001C TFC3 SGDID:S0000001, Chr I from 152168-146596, reverse complement, Verified ORF 
##                                                                                    5573 
##                     YAL002W VPS8 SGDID:S0000002, Chr I from 142709-148533, Verified ORF 
##                                                                                    5825 
##                     YAL003W EFB1 SGDID:S0000003, Chr I from 141176-144162, Verified ORF 
##                                                                                    2987 
## YAL005C SSA1 SGDID:S0000004, Chr I from 142433-138505, reverse complement, Verified ORF 
##                                                                                    3929 
## YAL007C ERP2 SGDID:S0000005, Chr I from 139347-136700, reverse complement, Verified ORF 
##                                                                                    2648 
##                    YAL008W FUN14 SGDID:S0000006, Chr I from 135916-138512, Verified ORF 
##                                                                                    2597 
##                     YAL009W SPO7 SGDID:S0000007, Chr I from 134856-137635, Verified ORF 
##                                                                                    2780
> x <- readDNAStringSet(filepath)
> x
##   A DNAStringSet instance of length 7
##     width seq                                          names               
## [1]  5573 ACTTGTAAATATATCTTTTAT...ATCGACCTTATTGTTGATAT YAL001C TFC3 SGDI...
## [2]  5825 TTCCAAGGCCGATGAATTCGA...AAATTTTTTTCTATTCTCTT YAL002W VPS8 SGDI...
## [3]  2987 CTTCATGTCAGCCTGCACTTC...TACTCATGTAGCTGCCTCAT YAL003W EFB1 SGDI...
## [4]  3929 CACTCATATCGGGGGTCTTAC...CCCGAAACACGAAAAAGTAC YAL005C SSA1 SGDI...
## [5]  2648 AGAGAAAGAGTTTCACTTCTT...TAATTTATGTGTGAACATAG YAL007C ERP2 SGDI...
## [6]  2597 GTGTCCGGGCCTCGCAGGCGT...TTTTGGCAGAATGTACTTTT YAL008W FUN14 SGD...
## [7]  2780 CAAGATAATGTCAAAGTTAGT...AAGGAAGAAAAAAAAATCAC YAL009W SPO7 SGDI...
> out1 <- tempfile()
> writeXStringSet(x, out1)
> ## 读写fastq文件
> filepath <- system.file("extdata", "s_1_sequence.txt", package = "Biostrings")
> fastq.geometry(filepath)
## [1] 256  36
> y <- readDNAStringSet(filepath, format = "fastq")
> out2 <- tempfile()
> writeXStringSet(y, out2, format = "fastq")
> ## bam文件
> library(Rsamtools)
> fl <- system.file("extdata", "ex1.bam", package = "Rsamtools", mustWork = TRUE)
> which <- RangesList(seq1 = IRanges(1000, 2000), seq2 = IRanges(c(100, 1000), 
+     c(1000, 2000)))
> p1 <- ScanBamParam(which = which)
> res1 <- scanBam(fl, param = p1)
## Warning: no BAM fields selected for input (niether 'bamWhat(param)' nor
##   'bamTag(param)' defined)
> names(res1)
## [1] "seq1:1000-2000" "seq2:100-1000"  "seq2:1000-2000"
> names(res1[[2]])
## character(0)
> ## GFF, BED, bedGraph, WIG, BigWig文件
> library(rtracklayer)
> track <- import(system.file("tests", "v1.gff", package = "rtracklayer"))

由此可见,Bioconductor对于生物数据的处理是非常轻松的。

2 thoughts on “生物信息学生R入门教程–输入与输出

  1. 您好:我想请教您怎么能从零基础学好R并应用于生物信息学中?我自己看了很久收获特别少

发表评论

电子邮件地址不会被公开。 必填项已用*标注