- fastp是一個(gè)非常方便的基因組數(shù)據(jù)處理軟件，之前在質(zhì)控的環(huán)節(jié)就有用到它。
- 這個(gè)軟件比較方便的地方就是對(duì)于處理前后的數(shù)據(jù)情況會(huì)輸出一個(gè)可視化的網(wǎng)頁(yè)報(bào)告，可以直接查看，報(bào)告內(nèi)容也十分的簡(jiǎn)單易懂。
- 對(duì)于單個(gè)數(shù)據(jù)或者少數(shù)數(shù)據(jù)直接這么查看是比較方便的，但是對(duì)于大批量的數(shù)據(jù)進(jìn)行處理以后，要還是一個(gè)一個(gè)報(bào)告查看，就讓人很抓狂了。
- 最近就碰到了這么一個(gè)問(wèn)題，于是寫(xiě)了一個(gè)python小腳本，進(jìn)行批量結(jié)果的可視化，做一下筆記，怕以后忘了。
- 要進(jìn)行批量可視化的前提是可視化的數(shù)據(jù)，fastp在處理完數(shù)據(jù)后，除了一個(gè)網(wǎng)頁(yè)格式報(bào)告之外，還有一個(gè)json格式的報(bào)告，這就是我們的數(shù)據(jù)來(lái)源。
- 本腳本使用python3.x的版本，設(shè)計(jì)第三方庫(kù)如下：
- pandas
- matplotlib
-
首先為這個(gè)腳本新建一個(gè)目錄，并在目錄下新建一個(gè)report目錄（目錄名隨意）和一個(gè)python腳本（腳本名隨意），report目錄用來(lái)放fastp的json格式報(bào)告，如圖，我的目錄下面就放了兩百多份報(bào)告。
目錄
- 這個(gè)實(shí)現(xiàn)本身很簡(jiǎn)單，主要要考慮的一點(diǎn)是WGS數(shù)據(jù)很多時(shí)候不是一個(gè)sample一對(duì)fastq文件，很多時(shí)候一個(gè)sample對(duì)應(yīng)好多對(duì)fastq文件，那么fastp也會(huì)輸出多個(gè)報(bào)告，對(duì)于這些報(bào)告我們要進(jìn)行合并。
- 下面就來(lái)解析一下腳本
先上完整腳本
import json, os
import pandas as pd
import matplotlib.pyplot as plt
result_list = os.listdir('./report')
result_dict = {}
merge_result_dict = {}
for i in result_list:
if i.endswith('.json'):
with open('./report/' + i, 'r') as f:
result_dict[i] = json.load(f)
for k,v in result_dict.items():
key = k.split('_')[0]
merge_result_dict[key] = {'total_bases': v['summary']['before_filtering']['total_bases'] + merge_result_dict.get(key, {'total_bases': 0}).get('total_bases', 0),
'q20_bases': v['summary']['before_filtering']['q20_bases'] + merge_result_dict.get(key, {'q20_bases': 0}).get('q20_bases', 0),
'q30_bases': v['summary']['before_filtering']['q30_bases'] + merge_result_dict.get(key, {'q30_bases': 0}).get('q30_bases', 0),
'total_reads': v['summary']['before_filtering']['total_reads']+ merge_result_dict.get(key, {'total_reads': 0}).get('total_reads', 0),
'dup': v['summary']['before_filtering']['total_reads'] * v['duplication']['rate'] + merge_result_dict.get(key, {'dup': 0}).get('dup', 0)}
df_ = pd.DataFrame(merge_result_dict).T
#原始數(shù)據(jù)量
plt.figure(1)
plt.scatter(df_.index, df_.total_bases)
plt.xlabel("sample name")
plt.ylabel("total_bases")
plt.title("Raw date: total_bases")
plt.plot(df_.index, [9e+10 for i in range(df_.shape[0])], 'r')
#q20
plt.figure(2)
plt.scatter(df_.index, df_.q20_bases / df_.total_bases)
plt.xlabel("sample name")
plt.ylabel("q20 percent")
plt.title("Raw date: q20 content")
#q30
plt.figure(3)
plt.scatter(df_.index, df_.q30_bases / df_.total_bases)
plt.xlabel("sample name")
plt.ylabel("q30 percent")
plt.title("Raw date: q30 content")
#dup
plt.figure(4)
plt.scatter(df_.index, df_.dup / df_.total_reads)
plt.xlabel("sample name")
plt.ylabel("duplication rate")
plt.title("Raw date: duplication rate")
#cut dup
plt.figure(5)
plt.scatter(df_.index, (1 - df_.dup / df_.total_reads) * df_.total_bases)
plt.xlabel("sample name")
plt.ylabel("bases")
plt.title("cut dup: bases")
plt.show()
-
首先導(dǎo)入相應(yīng)的庫(kù)，其中json和os是內(nèi)置庫(kù)，分別用來(lái)解析json文件以及查看文件。json格式的解析很簡(jiǎn)單，其實(shí)json可以看成是python里面的字典，下面是簡(jiǎn)單的json庫(kù)使用：
import json # 輸出 json 數(shù)據(jù) with open('data.json', 'w') as f: json.dump(data, f) # 讀取數(shù)據(jù) with open('data.json', 'r') as f: data = json.load(f)
一般情況下，在當(dāng)前目錄運(yùn)行的腳本，直接用相對(duì)路徑就行，但是也有可能你的報(bào)告不在當(dāng)前目錄下，那么第5行和第10行的
'./report'就要改成你報(bào)告所在目錄的位置，同時(shí)report這個(gè)目錄名也要改成你自己存放json格式的目錄名。-
如何合并一個(gè)sample的報(bào)告結(jié)果
for k,v in result_dict.items(): key = k.split('_')[0] ...這一步的工作就是合并結(jié)果，在我的結(jié)果里，雖然一個(gè)sample有時(shí)候會(huì)被分為好多個(gè)數(shù)據(jù)，但是數(shù)據(jù)的前綴的一部分都是一樣的，且以
_分割。 這個(gè)腳本提取了raw data中
total_bases和q30_bases的數(shù)據(jù)，這是因?yàn)檫@一塊是需要和測(cè)序公司核對(duì)的部分，一般公司需要交付 90G 的raw data，且Q30要80%以上。其實(shí)fastp結(jié)果里還有很多其他可以提取的數(shù)據(jù)，完全可以根據(jù)需要提取作圖，還是很簡(jiǎn)單的。最后一個(gè)圖是看一下，如果去掉duplicate的數(shù)據(jù)，那么raw data還有多少數(shù)據(jù)量。
展示一個(gè)結(jié)果圖。
-
腳本本身超級(jí)簡(jiǎn)單，我看半天實(shí)在也不知道還有哪里需要解釋的，如果有什么錯(cuò)誤或者寫(xiě)的不好的還請(qǐng)見(jiàn)諒了。
水平有限，要是存在什么錯(cuò)誤請(qǐng)指出，可發(fā)送郵箱至shiyuant@outlook.com！請(qǐng)大家多多批評(píng)指正，相互交流，共同成長(zhǎng)，謝謝?。?！
