Download Bioinformatics Toolbox User's Guide

2
High-Throughput Sequence Analysis
Load the human chromosome 9 from the reference file hs37.fasta. For this example, it
is assumed that you recovered the reference from the Bowtie indices using the bowtieinspect command; therefore hs37.fasta contains all the human chromosomes. To
load only the chromosome 9 you can use the option nave-value pair BLOCKREAD with the
fastaread function.
chr9 = fastaread('hs37.fasta','blockread',9)
chr9 =
Header: 'gi|224589821|ref|NC_000009.11| Homo sapiens chromosome 9, G...'
Sequence: 'NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN...'
Use the cpgisland function to find the CpG clusters. Using the standard definition for
CpG islands [4], 200 or more bp islands with 60% or greater CpGobserved/CpGexpected
ratio, leads to 1682 GpG islands found in chromosome 9.
cpgi = cpgisland(chr9.Sequence)
cpgi =
Starts: [1x1682 double]
Stops: [1x1682 double]
Use the getCounts method to calculate the ratio of aligned bases that are inside CpG
islands. For the first replicate of the sample HCT116, the ratio is close to 45%.
aligned_bases_in_CpG_islands = getCounts(bm_hct116_1,cpgi.Starts,cpgi.Stops,'method','s
aligned_bases_total = getCounts(bm_hct116_1,1,n,'method','sum')
ratio = aligned_bases_in_CpG_islands ./ aligned_bases_total
aligned_bases_in_CpG_islands =
1724363
aligned_bases_total =
3822804
2-86