经常有朋友问起自己要做什么分析，推荐一个电脑的配置。通常限制程序运行的最主要因素是内存，内存不足程序会直接运行不起来，CPU性能弱顶多是运行的慢，硬盘比较便宜，不需要特别评估。
针对这个问题，我们准备推出一系列测试推文，统计计算常用软件的运行时间、所需的最大物理内存（后面统计的都是所需最大物理内存，也就是配电脑所需的内存；程序运行过程中可能还会使用虚拟内存，这个参数也有统计，但不在分析范围之内）、生成文件所占的磁盘空间，以给后续选择做参考依据。
资源统计主要基于Time除了监控程序运行时间还能干这个？提到的time函数。
测试STAR构建基因组索引所需的计算资源
基因组序列和基因注释文件下载和整理
#下载当前最新版本的基因组wget-chttp://ftp.ensembl.org/pub/release-104/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz#处理基因组序列文件，染色体名字增加chr#且名字简化，去除无关信息gunzipHomo_sapiens.GRCh38.dna.primary_assembly.fa.gz-c|sed's/^>/>chr/'|cut-d''-f1>GRCh38.fa#下载当前最新基因注释文件wget-chttp://ftp.ensembl.org/pub/release-104/gtf/homo_sapiens/Homo_sapiens.GRCh38.104.gtf.gz#给染色体名字增加chrgunzipHomo_sapiens.GRCh38.104.gtf.gz-c|grep-v'^#'|sed's/^/chr/'>GRCh38.gtf
测试构建完整人类基因组所需的资源
调用STAR程序构建基因组索引
/usr/bin/time-oSTAR_human_genome.log-vSTAR--runModegenomeGenerate--runThreadN2--genomeDirstar_GRCh38--genomeFastaFilesGRCh38.fa--sjdbGTFfileGRCh38.gtf>STAR_human_genome_run.log2>&1
统计信息输出如下(这里只关注下面的3个统计信息，CPU利用率、运行耗时、最大物理内存)：
PercentofCPUthisjobgot:174%Elapsed(wallclock)time(h:mm:ssorm:ss):1:36:08Maximumresidentsetsize(kbytes):36945552
将人类基因组按染色体拆分模拟不同大小基因组构建索引的计算资源需求
采用染色体累加的方式，不断模拟不同大小的基因组对计算资源的需求。
/bin/rm-fGRCh38_tmp.seqidforiin`seq122`;do#累加染色体，第一次循环测试chr1，#第二次循环测试chr1，chr2echo"chr$i">>GRCh38_tmp.seqid#调用seqkit提取序列seqkitgrep-fGRCh38_tmp.seqidGRCh38.fa-oGRCh38_tmp.fa#提取对应的gtf注释awk'BEGIN{OFS=FS="\t"}ARGIND==1{id[$1]=1}ARGIND==2{if(id[$1]==1)print$0}'GRCh38_tmp.seqidGRCh38.gtf>GRCh38_tmp.gtf#统计构建基因组索引所需的计算资源/usr/bin/time-oSTAR_human_genome_${i}.log-vSTAR--runModegenomeGenerate--runThreadN2--genomeDirstar_GRCh38_${i}--genomeFastaFilesGRCh38_tmp.fa--sjdbGTFfileGRCh38_tmp.gtf>STAR_human_genome_${i}.run_log2>&1done统计输出结果同上，比如对chr1构成的基因组统计结果如下：
PercentofCPUthisjobgot:91%Elapsed(wallclock)time(h:mm:ssorm:ss):5:33.95Maximumresidentsetsize(kbytes):7719688
整合数据生成最终测试结果数据集
首先写一个awk脚本，整理并转换CPU使用率、程序耗时、最大物理内存需求。这个脚本存储为timeIntegrate.awk。
#!/usr/bin/awk-f#定义一个世界转换函数functionto_minutes(time_str){a=split(time_str,array1,":");minutes=0;count=1;for(i=a;i>=1;--i){minutes+=array1[count]*60^(i-2);count+=1;}returnminutes;}#设定特殊的输入分隔符，注意冒号后面的空格BEGIN{OFS="\t";FS=":";}#数据量计量ARGIND==1{datasize=$0;}#解析time的输出ARGIND==2{if(FNR==1&&outputHeader==1)print"Time_cost\tMemory_cost\tnCPU\tData_size";if($1~/Elapsed/){#时间转为分钟time_cost=to_minutes($2);}elseif($1~/Maximumresidentsetsize/){#内存转为Gmemory_cost=$2/10^6;}elseif($1~/CPU/){#cpu计算为n核cpu=($2+0)/100};}END{printtime_cost,memory_cost,cpu,datasize}使用上面的脚本，整合所有基因组数据构建所需的计算和时间资源信息。
/bin/rm-fGRCh38_star_genome_build.summaryforiin`seq122`;doawk-vi=${i}'BEGIN{for(j=1;j<=i;j++)scaffold["chr"j]=1;}\{if(scaffold[$1]==1)genomeSize+=$2}END{printgenomeSize/10^9}'star_GRCh38/chrNameLength.txt|\awk-voutputHeader=$i-f./timeIntegrate.awk-STAR_human_genome_${i}.log>>GRCh38_star_genome_build.summarydoneawk'{genomeSize+=$2}END{printgenomeSize/10^9}'star_GRCh38/chrNameLength.txt\|awk-f./timeIntegrate.awk-STAR_human_genome.log>>GRCh38_star_genome_build.summary最终获得结果文件如下：
Time_costMemory_costnCPUData_size5.565837.719690.910.24895611.612815.73510.930.4911514.754718.74951.090.68944615.904324.37671.340.8796623.360830.44861.641.061228.57835.45941.491.23233.418533.45711.411.3913537.900532.82321.341.5364939.834235.4091.381.6748839.538334.74461.51.8086840.045834.57461.61.9437742.469734.45931.612.0770446.180333.4321.542.1914153.134233.97071.542.2984556.111732.47661.512.4004458.941731.66881.482.4907859.954331.47971.512.5740461.132.7431.542.654416932.9211.642.7130369.033333.6921.652.7774775.934.17491.72.8241889.266734.70631.762.87596.133336.94561.743.09975
部分构建好的索引因为硬盘不足删除了，后期才想着应该统计下，故做了手动添加：
Time_cost(minutes)Memory_cost(G)nCPUData_size(G)Index_size(G)5.565837.719690.910.2489563.611.612815.73510.930.49115NA14.754718.74951.090.6894467.515.904324.37671.340.87966NA23.360830.44861.641.06121128.57835.45941.491.232NA33.418533.45711.411.391351437.900532.82321.341.53649NA39.834235.4091.381.674881639.538334.74461.51.80868NA40.045834.57461.61.943771942.469734.45931.612.07704NA46.180333.4321.542.191412153.134233.97071.542.29845NA56.111732.47661.512.400442358.941731.66881.482.49078NA59.954331.47971.512.574042461.132.7431.542.65441NA6932.9211.642.713032569.033333.6921.652.77747NA75.934.17491.72.824182689.266734.70631.762.875NA96.133336.94561.743.0997528
构建索引的时间随数据量的变化
构建基因组的索引所需时间跟基因组大小成正相关，大体分为3个阶段：
- 1.8G以内基因组构建索引的时间与基因组大小近乎完美的正相关。
- 1.8G-2.3G大小的基因组构建索引的时间基本相当，受基因组大小影响较小。
- 2.3G-3G大小的基因组构建索引的时间与基因组大小正相关，且时间随基因组大小变化的幅度大于基因组大小在1.8G以内时。
- 推测更大的基因组，时间需求也会更大。
此处只考虑了大小，但基因组序列组成的复杂性也会对构建索引有很大影响。
另外，后面应该尝试学习个模型做个简单预测。
#devtools::install_git("https://gitee.com/ct586/ImageGP")library(ImageGP)library(ggplot2)sp_scatterplot("GRCh38_star_genome_build.summary",melted=T,xvariable="Data_size",yvariable="Time_cost",smooth_method="auto",x_label="Genomesize(Gb)",y_label="Runningtime(minutes)")+scale_x_continuous(breaks=seq(0.5,3,by=0.5))+scale_y_continuous(breaks=seq(10,100,length=10))构建索引的所需内存随数据量的变化
- 在基因组大小超过1G时，STAR构建基因组索引需要的内存已经达到30G。
- 基因组大小在1.2G-2.8G时，STAR构建基因组索引需要的内存变化不大，在31G-35G之间徘徊，甚至出现了基因组变大，而所需内存变少的现象。
- （还不知道怎么解释这一现象）
- 基因组大小在3G时，需要的内存是37G。
- 推测导致内存增加的不只是因为数据量的增加，更因为数据变得更复杂了（这部分增加的基因组大小来源于性染色体和不能确定来源的基因组区域）
/usr/bin/time-oSTAR_human_genome.log-vSTAR--runModegenomeGenerate--runThreadN2--genomeDirstar_GRCh38--genomeFastaFilesGRCh38.fa--sjdbGTFfileGRCh38.gtf>STAR_human_genome_run.log2>&10
构建索引时对CPU的利用率
在构建索引时，我们只给程序分配了2个线程；检测结果也有些出乎意料，STAR似乎对多线程的利用并不充分，尤其时基因组比较小时，可能更多还是IO密集型的操作。
从趋势来看，随着基因组大小变大，CPU的双线程利用率也是逐渐增高的。这也就是说，有时核也不需要给太多，没那么大数据量时，给了多核不过是心里安慰。
/usr/bin/time-oSTAR_human_genome.log-vSTAR--runModegenomeGenerate--runThreadN2--genomeDirstar_GRCh38--genomeFastaFilesGRCh38.fa--sjdbGTFfileGRCh38.gtf>STAR_human_genome_run.log2>&11
生成的索引大小根基因组大小的关系（磁盘需求）
最开始统计时忽略了这个信息，本来time是可以输出Filesystemoutputs信息的，但不知道在这怎么会都是0，可能是操作系统的问题。
简单的根据最终生成的索引大小来做个分析和判断，一个近乎完美的正相关，基因组越大，产生的索引越大，大约是8-9倍基因组的大小。一个完整的人类基因构建的STAR索引需要28G的硬盘空间。
/usr/bin/time-oSTAR_human_genome.log-vSTAR--runModegenomeGenerate--runThreadN2--genomeDirstar_GRCh38--genomeFastaFilesGRCh38.fa--sjdbGTFfileGRCh38.gtf>STAR_human_genome_run.log2>&12