經常有朋友問起自己要做什么分析，推薦一個電腦的配置。通常限制程序運行的最主要因素是內存，內存不足程序會直接運行不起來，CPU性能弱頂多是運行的慢，硬盤比較便宜，不需要特別評估。

針對這個問題，我們準備推出一系列測試推文，統計計算常用軟件的運行時間、所需的最大物理內存 （后面統計的都是所需最大物理內存，也就是配電腦所需的內存；程序運行過程中可能還會使用虛擬內存，這個參數也有統計，但不在分析范圍之內）、生成文件所占的磁盤空間，以給后續選擇做參考依據。

資源統計主要基于Time除了監控程序運行時間還能干這個？提到的time函數。

測試 STAR 構建基因組索引所需的計算資源

基因組序列和基因注釋文件下載和整理

# 下載當前最新版本的基因組 wget -c http://ftp.ensembl.org/pub/release-104/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz # 處理基因組序列文件，染色體名字增加 chr # 且名字簡化，去除無關信息 gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz -c | sed 's/^>/>chr/' | cut -d ' ' -f 1 > GRCh38.fa # 下載當前最新基因注釋文件 wget -c http://ftp.ensembl.org/pub/release-104/gtf/homo_sapiens/Homo_sapiens.GRCh38.104.gtf.gz # 給染色體名字增加 chr gunzip Homo_sapiens.GRCh38.104.gtf.gz -c | grep -v '^#' | sed 's/^/chr/' >GRCh38.gtf

測試構建完整人類基因組所需的資源

調用 STAR 程序構建基因組索引

/usr/bin/time -o STAR_human_genome.log -v STAR --runMode genomeGenerate --runThreadN 2 --genomeDir star_GRCh38 --genomeFastaFiles GRCh38.fa --sjdbGTFfile GRCh38.gtf > STAR_human_genome_run.log 2>&1

統計信息輸出如下(這里只關注下面的 3 個統計信息，CPU利用率、運行耗時、最大物理內存)：

Percent of CPU this job got: 174% Elapsed (wall clock) time (h:mm:ss or m:ss): 1:36:08 Maximum resident set size (kbytes): 36945552

將人類基因組按染色體拆分模擬不同大小基因組構建索引的計算資源需求

采用染色體累加的方式，不斷模擬不同大小的基因組對計算資源的需求。

/bin/rm -f GRCh38_tmp.seqid for i in `seq 1 22`; do# 累加染色體，第一次循環測試 chr1，# 第二次循環測試 chr1，chr2echo "chr$i" >>GRCh38_tmp.seqid# 調用 seqkit 提取序列seqkit grep -f GRCh38_tmp.seqid GRCh38.fa -o GRCh38_tmp.fa# 提取對應的 gtf 注釋awk 'BEGIN{OFS=FS="\t"}ARGIND==1{id[$1]=1}ARGIND==2{if(id[$1]==1) print $0}' GRCh38_tmp.seqid GRCh38.gtf >GRCh38_tmp.gtf# 統計構建基因組索引所需的計算資源/usr/bin/time -o STAR_human_genome_${i}.log -v STAR --runMode genomeGenerate --runThreadN 2 --genomeDir star_GRCh38_${i} --genomeFastaFiles GRCh38_tmp.fa --sjdbGTFfile GRCh38_tmp.gtf > STAR_human_genome_${i}.run_log 2>&1 done

統計輸出結果同上，比如對chr1構成的基因組統計結果如下：

Percent of CPU this job got: 91% Elapsed (wall clock) time (h:mm:ss or m:ss): 5:33.95 Maximum resident set size (kbytes): 7719688

整合數據生成最終測試結果數據集

首先寫一個awk腳本，整理并轉換 CPU 使用率、程序耗時、最大物理內存需求。這個腳本存儲為timeIntegrate.awk。

#!/usr/bin/awk -f# 定義一個世界轉換函數 function to_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;}return minutes;} # 設定特殊的輸入分隔符，注意冒號后面的空格 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);} else if($1~/Maximum resident set size/){# 內存轉為 Gmemory_cost=$2/10^6;} else if($1~/CPU/) {# cpu 計算為 n核cpu=($2+0)/100}; }END{print time_cost, memory_cost, cpu, datasize}

使用上面的腳本，整合所有基因組數據構建所需的計算和時間資源信息。

/bin/rm -f GRCh38_star_genome_build.summary for i in `seq 1 22`; doawk -v i=${i} 'BEGIN{for(j=1;j<=i;j++) scaffold["chr"j]=1;}\{if(scaffold[$1]==1) genomeSize+=$2}END{print genomeSize/10^9}' star_GRCh38/chrNameLength.txt | \awk -v outputHeader=$i -f ./timeIntegrate.awk - STAR_human_genome_${i}.log >>GRCh38_star_genome_build.summary doneawk '{genomeSize+=$2}END{print genomeSize/10^9}' star_GRCh38/chrNameLength.txt \ | awk -f ./timeIntegrate.awk - STAR_human_genome.log >>GRCh38_star_genome_build.summary

最終獲得結果文件如下：

Time_cost Memory_cost nCPU Data_size 5.56583 7.71969 0.91 0.248956 11.6128 15.7351 0.93 0.49115 14.7547 18.7495 1.09 0.689446 15.9043 24.3767 1.34 0.87966 23.3608 30.4486 1.64 1.0612 28.578 35.4594 1.49 1.232 33.4185 33.4571 1.41 1.39135 37.9005 32.8232 1.34 1.53649 39.8342 35.409 1.38 1.67488 39.5383 34.7446 1.5 1.80868 40.0458 34.5746 1.6 1.94377 42.4697 34.4593 1.61 2.07704 46.1803 33.432 1.54 2.19141 53.1342 33.9707 1.54 2.29845 56.1117 32.4766 1.51 2.40044 58.9417 31.6688 1.48 2.49078 59.9543 31.4797 1.51 2.57404 61.1 32.743 1.54 2.65441 69 32.921 1.64 2.71303 69.0333 33.692 1.65 2.77747 75.9 34.1749 1.7 2.82418 89.2667 34.7063 1.76 2.875 96.1333 36.9456 1.74 3.09975

部分構建好的索引因為硬盤不足刪除了，后期才想著應該統計下，故做了手動添加：

Time_cost (minutes) Memory_cost (G) nCPU Data_size (G) Index_size (G) 5.56583 7.71969 0.91 0.248956 3.6 11.6128 15.7351 0.93 0.49115 NA 14.7547 18.7495 1.09 0.689446 7.5 15.9043 24.3767 1.34 0.87966 NA 23.3608 30.4486 1.64 1.0612 11 28.578 35.4594 1.49 1.232 NA 33.4185 33.4571 1.41 1.39135 14 37.9005 32.8232 1.34 1.53649 NA 39.8342 35.409 1.38 1.67488 16 39.5383 34.7446 1.5 1.80868 NA 40.0458 34.5746 1.6 1.94377 19 42.4697 34.4593 1.61 2.07704 NA 46.1803 33.432 1.54 2.19141 21 53.1342 33.9707 1.54 2.29845 NA 56.1117 32.4766 1.51 2.40044 23 58.9417 31.6688 1.48 2.49078 NA 59.9543 31.4797 1.51 2.57404 24 61.1 32.743 1.54 2.65441 NA 69 32.921 1.64 2.71303 25 69.0333 33.692 1.65 2.77747 NA 75.9 34.1749 1.7 2.82418 26 89.2667 34.7063 1.76 2.875 NA 96.1333 36.9456 1.74 3.09975 28

構建索引的時間隨數據量的變化

構建基因組的索引所需時間跟基因組大小成正相關，大體分為 3 個階段：

1.8 G以內基因組構建索引的時間與基因組大小近乎完美的正相關。

1.8 G - 2.3 G大小的基因組構建索引的時間基本相當，受基因組大小影響較小。

2.3 G - 3 G大小的基因組構建索引的時間與基因組大小正相關，且時間隨基因組大小變化的幅度大于基因組大小在1.8 G以內時。

推測更大的基因組，時間需求也會更大。

此處只考慮了大小，但基因組序列組成的復雜性也會對構建索引有很大影響。

另外，后面應該嘗試學習個模型做個簡單預測。

# 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 ="Genome size (Gb)", y_label = "Running time (minutes)") +scale_x_continuous(breaks=seq(0.5,3, by=0.5)) +scale_y_continuous(breaks=seq(10,100,length=10))

構建索引的所需內存隨數據量的變化

在基因組大小超過1 G時，STAR 構建基因組索引需要的內存已經達到30 G。

基因組大小在1.2 G - 2.8 G時，STAR 構建基因組索引需要的內存變化不大，在31 G-35 G之間徘徊，甚至出現了基因組變大，而所需內存變少的現象。

（還不知道怎么解釋這一現象）

基因組大小在3 G時，需要的內存是37 G。

推測導致內存增加的不只是因為數據量的增加，更因為數據變得更復雜了（這部分增加的基因組大小來源于性染色體和不能確定來源的基因組區域）

sp_scatterplot("GRCh38_star_genome_build.summary", melted = T, xvariable = "Data_size",yvariable = "Memory_cost", smooth_method = "auto",x_label ="Genome size (Gb)", y_label = "Maximum physical memory required (Gb)") +scale_x_continuous(breaks=seq(0.5,3, by=0.5)) +scale_y_continuous(breaks=seq(5,40, by=5))

構建索引時對 CPU 的利用率

在構建索引時，我們只給程序分配了 2 個線程；檢測結果也有些出乎意料，STAR 似乎對多線程的利用并不充分，尤其時基因組比較小時，可能更多還是 IO 密集型的操作。

從趨勢來看，隨著基因組大小變大，CPU 的雙線程利用率也是逐漸增高的。這也就是說，有時核也不需要給太多，沒那么大數據量時，給了多核不過是心里安慰。

sp_scatterplot("GRCh38_star_genome_build.summary", melted = T, xvariable = "Data_size",yvariable = "nCPU", smooth_method = "auto",x_label ="Genome size (Gb)", y_label = "Number of CPUs used") +scale_x_continuous(breaks=seq(0.5,3, by=0.5)) +scale_y_continuous(breaks=seq(0,2, by=0.1))

生成的索引大小根基因組大小的關系 （磁盤需求）

最開始統計時忽略了這個信息，本來time是可以輸出File system outputs信息的，但不知道在這怎么會都是0，可能是操作系統的問題。

簡單的根據最終生成的索引大小來做個分析和判斷，一個近乎完美的正相關，基因組越大，產生的索引越大，大約是8-9倍基因組的大小。一個完整的人類基因構建的 STAR 索引需要28 G的硬盤空間。

sp_scatterplot("GRCh38_star_genome_build.summary", melted = T, xvariable = "Data_size",yvariable = "Index_size", smooth_method = "auto",x_label ="Genome size (Gb)", y_label = "Disk space usages (Gb)") +scale_x_continuous(breaks=seq(0.5,3, by=0.5)) +scale_y_continuous(breaks=seq(2,28, by=2))

未完待續~~~~

后面就是不同數據量比對時所需的計算資源評估了！！！

AlphaFold2開源了，不是土豪也不會編程的你怎么蹭一波？

50T內存？百萬機時？頭一次見這么耗費內存和機時的分析？

往期精品(點擊圖片直達文字對應教程)

機器學習

后臺回復“生信寶典福利第一波”或點擊閱讀原文獲取教程合集

總結

以上是生活随笔為你收集整理的什么配置的电脑可满足基因组索引构建的需求？的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。