如果您想要快速又髒的東西，可以使用 samtools faidx 快速索引FASTA，然後在命令行上將R的lengths列放入R（其他語言可用）。 pre> samtools faidx $ fastacut -f2 $ fasta.fai | Rscript -e'data <- as.numeric（readLines（“ stdin”））;摘要（數據）； hist（data）'

這將輸出統計摘要，並在當前目錄Rplots.pdf中創建一個包含直方圖的PDF。

納米孔讀取的統計數據非常棘手，因為一次運行中可能會出現大量讀取長度。我發現最好的顯示長度的方法是在x軸（長度）和y軸上使用對​​數刻度（順序基數或計數，具體取決於首選項）。

我自己執行此操作的腳本：一個用於生成讀取的長度，另一個用於以各種方式繪製長度分佈。生成讀取長度的腳本還會將基本長度摘要統計信息吐出為標準錯誤：

  $〜/ scripts / fastx-length.pl > lengths_mtDNA_drawn.txt總序列：2110總長度：5.106649 Mb最長序列： 107.414 kb最短序列：219 b平均長度：2.42 kb中位數長度：1.504 kb N50：336個序列; L50：3.644 kb N90：1359個序列; L90：1.103 kb $〜/ scripts / length_plot.r lengths_mtDNA_drawn.txtlengths_mtDNA_called.txt ...完成序列號：2110長度分位數：0％10％20％30％40％50％60％70％219.0 506.9 724.4 953.0 1196.2 1503.0 1859.2 2347.3 80％90％100％3128.2 4804.7 107414.0  

以下是幾個生成的圖形：

可以在此處找到生成這些腳本的腳本：

• https：// github.com/gringer/bioinfscripts/blob/master/fastx-length.pl
• https://github.com/gringer/bioinfscripts/blob/master/length_plot.r

使用Biopython和matplotlib似乎確實可行，它實際上可以歸結為三行代碼來獲取該圖：

  import Bio，pandaslengths = map（len ，Bio.SeqIO.parse（'/ path / to / the / seqs.fasta'，'fasta'））pandas.Series（lengths）.hist（color ='gray'，bins = 1000） 

當然，您可能想要製作一個更長的腳本，該腳本可以從命令行調用，並帶有幾個選項。歡迎您使用我的：

 ＃！/ usr / bin / env python2“”“一個定制腳本來繪製長度在fasta文件中的分佈。由Lucas Sinclair.Kopimi撰寫。您可以從shell這樣使用該腳本：$ ./fastq_length_hist --input seqs.fasta --out seqs.pdf“”“ ##################### ################################################ ##########模塊＃從argparse導入argparse，sys，time，getpass，locale從Bio導入SeqIOimport熊貓＃Matplotlib #import matplotlibmatplotlib.use（'Agg'，warn = False）from matplotlib import pyplot＃ ################################################ ############################## desc =“ fasta_length_hist v1.0” parser = argparse.ArgumentParser（description = desc，formatter_class = RawTextHelpFormatter ）＃所有必需的參數＃parser.add_argument（“-input”，help =“要處理的fasta文件”，type = str）parser.add_argument（“-out”，type = str）＃所有可選參數＃parser.add_argument（“-x_log”，默認=真實，典型e = bool）parser.add_argument（“-y_log”，默認= True，類型= bool）＃解析#args = parser.parse_args（）input_path = args.inputoutput_path = args.outx_log = bool（args.x_log）y_log = bool（args.y_log）########################################### ####################################讀取#lengths = map（len，SeqIO.parse（ input_path，'fasta'））＃報告＃sys.stderr.write（“讀取所有長度（％i序列）\ n”％len（長度））sys.stderr.write（“最長序列：％i bp \ n” ％max（lengths））sys.stderr.write（“最短序列：％i bp \ n”％min（lengths））
sys.stderr.write（“製作圖形... \ n”）＃數據#values = pandas.Series（lengths）＃Plot #fig = pyplot.figure（）axes = values.hist（color ='gray'，bins = 1000）fig = pyplot.gcf（）title ='序列長度的分佈'axes.set_title（title）axes.set_xlabel（'序列中的核苷酸數'）axes.set_ylabel（'具有此長度的序列數'）axes .xaxis.grid（False）＃記錄#if x_log：axes.set_yscale（'symlog'）if y_log：axes.set_xscale（'symlog'）＃調整＃width = 18.0;身高= 10.0;底部= 0.1;最高= 0.93;左= 0.07; right = 0.98fig.set_figwidth（width）fig.set_figheight（height）fig.subplots_adjust（hspace = 0.0，bottom = bottom，top = top，left = left，right = right）＃數據和源＃fig.text（0.99， 0.98，time.asctime（），horizo​​ntalalignment ='right'）fig.text（0.01，0.98，'user：'+ getpass.getuser（），horizo​​ntalalignment ='left'）＃好的數字分組#sep =（'x' ，'y'）如果Sep中的'x'：locale.setlocale（locale.LC_ALL，``）seperate = lambda x，pos：locale.format（“％d”，x，grouping = True）axes.xaxis.set_major_formatter （matplotlib.ticker.FuncFormatter（seperate））如果Sep中的'y'：locale.setlocale（locale.LC_ALL，''）seperate = lambda x，pos：locale.format（“％d”，x，grouping = True） axes.yaxis.set_major_formatter（matplotlib.ticker.FuncFormatter（seperate））＃保存＃fig.savefig（output_path，format ='pdf'） 

EDIT-示例輸出：

有幾種潛在的方法。例如：Biopython教程中的

• 序列長度直方圖（來自基於Ruby的生物件框架）
• plot_distribution
• 獲得序列長度的各種解決方案，包括bioawk和 EMBOSS infoseq

其中哪些是“快速高效”使用10 GB的文件...這很難事先說清楚。您可能需要嘗試對其中一些進行基準測試。

bioawk對於這種任務可能相當有效。

  $ bioawk -c fastx'{histo [length（$ seq）] ++} END {對於（歷史上的l）打印l，histo [l]}'\ |排序-n0 332701 15422 11323 33974 87765 118846 124747 143418 131659 1546710 2108911 3046912 4520413 6231114 8874415 11576716 14077017 19181018 31308819 51811120 109786721 472971522 657555723 273406224 101547625 49332326 32382727 16441928 10712029 134938 23674637 24938234 63824 847238 63824 847238 18423 847 23184-2324 23268-234 23264-234 23268-234 63824-234-387-184 23加大碼184拍打瀏覽器5948 5049 5250 4851 4252 3953 2854 4955 5956 5557 5158 5559 4360 5261 5662 4863 6764 9565 488  

-c fastx 告訴程序將數據解析為fastq或Fasta。這樣就可以訪問記錄的不同部分，如 $ name ， $ seq （對於fastq格式為 $ qual ）。 awk代碼（bioawk基於awk，因此您可以使用 awk中想要的任何語言功能）。

在單引號之間包含一系列 <condition> {<action> } 塊。

第一個沒有 <condition> 部分，這意味著將對每個記錄執行該部分。在這裡，它更新了我命名為“ histo”的表中的長度計數。 length 是awk中的預定義函數。

在第二個塊中， END 條件意味著我們希望在所有輸入都已完成之後執行它處理。動作部分包括循環記錄的長度值，並將它們與關聯的計數一起打印。

將輸出通過管道傳遞到 sort -n 以便對結果進行數字排序。

在我的工作站上，上面的代碼花了20秒才能執行1.2G的fasta文件。

您曾特別詢問過FASTA文件，但是在評估高錯誤長時間讀取的數據時，務必共同考慮讀取的長度和質量，這一點很重要。 FASTA文件不提供質量。這些信息將幫助您確定運行是否成功，“高質量”的讀數有多少等。

我最初在此處發布了完整的答案，建議 pauvre，但由於您似乎只有FASTA文件，所以我認為這有點偏離主題。

我建議生成FASTQ文件，但不確定是否具有原始庫-稱為fast5文件。如果是這樣，請使用 poretools 如下生成FASTQ（ poretools文檔以生成FASTQ文件）：

  poretools fastq fast5 /  

然後我建議使用 pauvre讀取長度和質量，如下所示生成熱圖和直方圖裕度圖。

[注意：我不是pauvre的原始作者，但我現在正在為這個項目做貢獻]

這並不是您所要的，但是您可以使用 poretools