使用全基因组数据
这个例子展示了如何创建一个内存映射文件序列数据和使用它没有所有的基因组序列加载到内存中。整个基因组可用于人类,老鼠,老鼠,河豚和其他生物模型。对于许多这些生物体的染色体可以几亿碱基对长。处理如此庞大的数据集可以挑战你可能遇到的限制您使用的硬件和软件。这个例子展示了一个方法解决这些限制在MATLAB®。
大数据集的处理问题
解决技术计算问题,需要处理和分析大量数据的计算机系统上提出了高要求。大型数据集占用大量的内存中处理和需要许多操作来计算一个解决方案。它也需要很长时间才能从大型数据文件访问信息。
然而,计算机系统有限内存和CPU速度有限。可用资源不同的处理器和操作系统,后者也消耗资源。例如:
32位处理器和操作系统可以解决2 ^ 32 = 4294967296 = 4 GB的内存(也称为虚拟地址空间)。Windows®XP和Windows®2000只有2 GB的虚拟内存分配给每个进程(如MATLAB)。在UNIX®上,虚拟内存分配给过程system-configurable,通常在3 GB。MATLAB应用程序进行计算,例如,可能需要存储除了用户的任务。主要的问题在处理大量数据是程序的内存需求可以超过可用的平台。例如,MATLAB生成一个“内存溢出”错误当数据需求超过约1.7 GB在Windows XP。
为更多的细节在内存管理和大型数据集,明白了性能和内存。
在典型的32位机器上,一个数据集的最大大小,您可以使用MATLAB几百MB,或大小的大型染色体。内存映射文件允许MATLAB解决这个限制并使您能够使用非常大的数据集在一个直观的方式。
全基因组数据集
最近的全基因组数据集可以从下载运用网站。提供的数据是在几个格式。这些都是定期更新新序列信息。这个例子将使用人类DNA FASTA格式存储的数据。1号染色体(GRCh37.56发布的2009年9月)65.6 MB的压缩文件。解压缩后的文件大约250 mb。MATLAB使用每字符2字节,如果文件读入MATLAB,它需要大约500 mb的内存。
这个示例假设您已经下载并解压FASTA文件到你的本地目录。改变变量的名称FASTAfilename如果合适。
FASTAfilename =“Homo_sapiens.GRCh37.56.dna.chromosome.1.fa”;fileInfo = dir ((FASTAfilename))
fileInfo =结构体字段:名称:“Homo_sapiens.GRCh37.56.dna.chromosome.1。fa的文件夹:/ mathworks /中心/ qe / test_data / Bioinformatics_Toolbox / v000 / demoData biomemorymapdemo的日期:01 - 2月2013 11:54:41字节:253404851 isdir: 0 datenum: 7.3527 e + 05
内存映射文件
内存映射允许MATLAB访问数据文件好像是在内存中。您可以使用标准的MATLAB索引操作来访问数据。看到的文档memmapfile为更多的细节。
你可以映射FASTA文件和直接访问数据。然而FASTA格式文件包括新行字符。的memmapfile函数以同样的方式对待这些字符作为所有其他字符。删除这些内存映射文件之前将索引操作简单。同时,内存映射不直接处理字符数据,所以你要把数据作为8位整数(uint8类)。这个函数nt2int在生物信息学工具箱™可用于字符信息转换为整数值。int2nt用于转换回字符。
首先打开FASTA文件和提取标题。
fidIn = fopen (FASTAfilename,“r”);头= fgetl (fidIn)
头= ' > 1 dna:染色体染色体:GRCh37:1:1:249250621:1”
打开文件内存映射。
[fullPath,文件名,扩展]= fileparts (FASTAfilename);mmFilename =[文件名“功能”]fidOut = fopen (mmFilename,' w ');
mmFilename = ' Homo_sapiens.GRCh37.56.dna.chromosome.1.mm '
阅读FASTA 1 mb的文件,删除新行字符,转换uint8,写信给MM文件。
换行符= sprintf (' \ n ');blockSize = 2 ^ 20;而~ feof (fidIn)%读入数据charData =从文件中读(fidIn blockSize,“*字符”)”;%去除新行charData = strrep (charData换行,”);%转换为整数intData = nt2int (charData);%写入新文件写入文件(fidOut intData,“uint8”);结束
关闭文件。
文件关闭(fidIn);文件关闭(fidOut);
新文件大小相同的旧文件但不包含新行或消息头信息。
mmfileInfo = dir (mmFilename)
mmfileInfo =结构体字段:名称:“Homo_sapiens.GRCh37.56.dna.chromosome.1。毫米的文件夹:/ tmp / Bdoc22a_1891349_66459 / tp4e06d876 ex57563178的日期:2月26日—- 2022 12:38:14字节:249250621 isdir: 0 datenum: 7.3858 e + 05
访问内存映射文件中的数据
命令memmapfile构造一个memmapfile对象,将新文件映射到内存中。为了做到这一点,它需要知道文件的格式。这个文件的格式很简单,但更复杂的格式可以映射。
chr1 = memmapfile (mmFilename,“格式”,“uint8”)
chr1 =文件名:“/ tmp / Bdoc22a_1891349_66459 / tp4e06d876 ex57563178 / Homo_sapiens.GRCh37.56.dna.chromosome.1。毫米的可写:假抵消:0格式:uint8的重复:正数据:249250621 x1 uint8数组
MEMMAPFILE对象
memmapfile对象的各种属性。文件名存储文件的完整路径。可写的显示的数据是否可以修改。请注意,如果你修改数据,这也将修改原始文件。抵消允许您指定使用的空间任何头信息。格式表示数据格式。重复用于指定有多少块(所定义的格式)地图。这可能是有用的限制多少内存是用来创建内存映射。这些属性可以相同的方式访问其他MATLAB数据。更多细节见类型帮助memmapfile或医生memmapfile。
chr1.Data (1:10)
ans = 10 x1 uint8列向量15个15个15个15个15个15个15个15个15个15
您可以使用索引访问任何地区的数据操作。
chr1.Data (10000000:10000010) '
ans = 1 x11 uint8行向量1 1 2 2 2 2 3 4 2 4 2
记住,核苷酸信息转换为整数。您可以使用int2nt拿回的序列信息。
int2nt (chr1.Data (10000000:10000010)”)
ans = ' AACCCCGTCTC '
或使用seqdisp显示顺序。
seqdisp (chr1.Data (10000000:10001000)”)
ans = 17 x71 char数组' 1 AACCCCGTCT CTACAATAAA TTAAAATATT AGCTGGGCAT GGTGGTGTGT GCTTGTAGTC”“61年CCAGCTACTT GGCGGGCTGA GGTGGGAGAA TCATCCAAGC CTTGGAGGCA GAGGTTGCAG”“121年TGAGCTGAGA TTGTGACACT GCACTCCAGC CTGGGAGACA GAGTGAGACT CCTACTCAAA”“181年AAAAAACAAA AAACAAAAAA CAAACCACAA AACTTTCCAG GTAACTTATT AAAACATGTT”“241年TTTTGTTTGT TTTGAGACAG AGTCTTGCTC TGTCGCCCAG GCTGGAGTGC AGTGGAGCAA”“301年TCTCAGCTCA CTGCAAGCTC CGCCTCCCGG GTTCACACCA TTCTCCTGCC TCAGCCTCCC”“361年GAGTAGCTAG GACTATAGGC ACCCGCCACC ACGCCCAGCT TATTTTTTTT GTATTTTTTA”“421年GTAGAGACGG GGTTTCATCG TGTTAGCCAG GATGGTCTCG ATCTCCTGAC CTCGTGATCC”“481年GCCCACCTCA GCCTCCCAAA GTGCTGGGAT TACAGGCGTG AGCCACTGCA CCCGGCCTAG”“541年TTTTTGTATA TTTTTTTTAG TAGAGACAGG GTTTCACCAT GTTAGCCAGG ATGGTCTCAA”“601年TCTCCTGACC TCGTGATCCG CCCGCCTCGG CCTCCCAAAG TGCTGGGGTT ACAGGCGTGA”“661年GCCACCGCAC ACAGCATTAA AGCATGTTTT ATTTTCCTAC ACATAATGAA ATCATTACCA”“721年GATGATTTGA CATGTGTACT TCATTGGAGA GGATTCTTAC AGTATATTCA AAATTAAATA”“781年TAATGACAAA AAATTACTAC CTAATCTATT AAAATTGGCA TAAGTCATCT ATGATCATTA”“841年ATGATATGCA AACATAAACA AGTATTATAC CCAGAAGTGT AATTTATTGT AGCTACATCT”“901年TATGTATAAT AGTTTAGTGG ATTTTTCCTG GAAATTGTCC ATTTTAATTT TTCTCTTAAG”“961年TCTGTGGAAT TTTCCAGTAA AAGTCAAGGC AAACCCAAGA T '
分析整个染色体
现在,您可以很容易地访问整个染色体,可以分析数据。这个例子展示了一种方法看GC含量沿染色体。
您提取的500000块bp和计算GC含量。
计算有多少块使用。
numNT =元素个数(chr1.Data);blockSize = 500000;numBlocks =地板(numNT / blockSize);
一种帮助MATLAB性能在处理大型数据集时是“preallocate”空间数据。这允许MATLAB来分配足够的空间的所有数据,而不必在小块生长数组。这将加快速度,也保护你的问题太大存储的数据。预分配数组的更多细节,请参阅://www.tianjin-qmedu.com/matlabcentral/answers/99124-how-do-i-pre-allocate-memory-when-using-matlab
一个简单的方法来preallocate是使用数组0函数。
率= 0 (numBlocks + 1, - 1);
遍历寻找C和G的数据,然后把这个数字除以总数,T、C、G .这将需要大约一分钟。
一个= nt2int (“一个”);C = nt2int (“C”);G = nt2int (‘G’);T = nt2int (“T”);为数= 1:numBlocks%计算块的指标开始= 1 + blockSize *(把1);停止= blockSize *计数;%提取块块= chr1.Data(启动,停止);%找到GC和内容gc =(总和(块= = G |块= = C));=(总和(块= = |块= = T));%计算GC总比已知的核苷酸比(数)= gc / (gc +);结束
最后一块小所以把这作为一种特殊情况。
块= chr1.Data(停止+ 1:结束);gc =(总和(块= = G |块= = C));=(总和(块= = |块= = T));比(结束)= gc / (gc +);
GC的情节智人染色体1的内容
xAxis = [1: blockSize: numBlocks * blockSize numNT);情节(xAxis比率)包含(“碱基对”);ylabel (的相对GC内容”);标题(的相对GC含量智人染色体1 ')
该地区的中心情节大约140 mbp Ns是一个大型的地区。
seqdisp (chr1.Data (140000000:140001000))
ans = 17 x71 char数组' 1 NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN”“61年NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN”“121年NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN”“181年NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN”“241年NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN”“301年NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN”“361年NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN”“421年NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN”“481年NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN”“541年NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN”“601年NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN”“661年NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN”“721年NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN”“781年NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN”“841年NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN”“901年NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN”“961年NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN NNNNNNNNNN N '
寻找GC含量高的地区
您可以使用找到确定GC含量高的地区。
指数=找到(比> 0.5);范围= [(1 + blockSize *(索引设置为1)),blockSize *指数);流(“地区% d: % d GC含量% f \ n ',范围,比率(指数)”)
区域500001:1000000 GC含量0.501412区域1000001:1500000 GC含量0.598332区域1500001:2000000 GC含量0.539498区域2000001:2500000 GC含量0.594508区域2500001:3000000 GC含量0.568620区域3000001:3500000 GC含量0.584572区域3500001:4000000 GC含量0.548137区域6000001:6500000 GC含量0.545072区域9000001:9500000 GC含量0.506692区域10500001:11000000 GC含量0.511386区域11500001:12000000 GC含量0.519874区域16000001:16500000 GC含量0.513082区域17500001:18000000 GC含量0.513392区域21500001:22000000 GC含量0.505598区域156000001:156500000 GC含量0.504446区域156500001:157000000 GC含量0.504090区域201000001:201500000 GC含量0.502976区域228000001:228500000 GC含量0.511946
如果你想删除临时文件,您必须首先清除memmapfile对象。
清晰的chr1删除(mmFilename)
