miRNA是一類長約22 nt的非編碼小RNA。miRNA通過堿基互補配對結(jié)合到靶mRNAs 3’UTRs以調(diào)控數(shù)百個基因的表達。因此，miRNA參與幾乎所有生命活動，與生物體的生長、發(fā)育、疾病、衰老和死亡息息相關(guān)，可以作為疾病診斷和檢測的生物標志物。

日常數(shù)據(jù)分析中，我們經(jīng)常會遇到各種來源的miRNA測序原始數(shù)據(jù)，有些測36 bp，有些測50 bp，有些測75 bp或者更長的。形形色色的reads長度，adaptor接頭序列，眾多測序公司來源的數(shù)據(jù)，給數(shù)據(jù)分析帶來了諸多不便。其中比較常見的問題就是：接頭序列是什么？

小編就曾經(jīng)遇到過這個問題，問某公司：“你們測序的接頭序列是什么”？答復(fù)：“保密”！小編當時就想懟他：不知道接頭序列，就沒法把接頭去掉，不去掉接頭怎么往下分析？都是一個圈子混的，保不保密我不知道嗎？

今天小編就教你如何在不知道接頭序列的情況下，把接頭序列從原始fastq文件中扒出來，同時把read的組成結(jié)構(gòu)也分析一下。

1，fastq數(shù)據(jù)格式

圖1. fastq格式

Fastq格式是測序原始數(shù)據(jù)存儲標準格式，存儲了每條read的序列及堿基質(zhì)量分數(shù)。如圖1所示，每4行代表一條read。

以Read1為例：

第1行代表序列的編號，以@開頭；一般包括：機器編碼（A00212），run id（123），flowcell編號（H3TTTDSX3），lane號（4），flowcell上的tile編號（1101），X坐標（8919）和Y坐標（1000），其后是barcode信息，用于將多個樣品拆開（demultiplex）。

第2行是read序列，一般ATCG，有時候有N。N表示熒光信號干擾無法判斷到底是什么堿基。

第3行是識別碼，以+開頭；后邊接序列標示符、描述信息等，或者干脆就是+

第4行是堿基質(zhì)量分數(shù)，與read序列一一對應(yīng)，即每個堿基有一個分數(shù)，表明該堿基的可信程度。由于不同型號機器出來的值不太一樣，因此可以根據(jù)這些值大概判斷測序儀型號。

2，miRNA read構(gòu)成

通常，miRNA read有三種構(gòu)成方式（如圖2所示）：

A：miRNA序列后邊直接連接adaptor接頭序列

B：miRNA序列前面加幾個堿基，后邊接miRNA序列，然后連接adaptor接頭序列

C：miRNA序列后邊加幾個堿基，然后連接adaptor接頭序列

圖2. miRNA read構(gòu)成

3，miRNA常見接頭序列

通常illumina測序的小RNA接頭序列有：TGGAATTCTCGGGTGCCAAGGAACTCC、AACTGTAGGCACCATCAAT、AGATCGGAAGAGCACACGTCT等，拿到原始數(shù)據(jù)后，首先可以先用這3種序列去搜下。搜到了就可以先用這些接頭去預(yù)分析下，拿到去接頭后的miRNA長度分布圖，如果長度分布正常（在22nt）附近有峰，則表示接頭序列正常去掉，這種一般是圖2中的A情況。對于B和C就無能為力了。

4，Clustalw檢測接頭序列

4.1 fastq collapse

在完全未知的情況下，我們可以這樣操作。首先，將fastq格式轉(zhuǎn)成fa格式，并且對每條不同的read序列進行計數(shù)，所謂的collapse步驟。

可以使用Galaxy在線工具進行操作。https://usegalaxy.org/，搜索collapse，上傳fa或者fastq文件即可，其功能是統(tǒng)計每條unique序列的個數(shù)

圖3. read序列collapse

以小編手里的fastq文件為例，經(jīng)過collapse后，前10條序列為：

圖4.top 10 read

肉眼可見，中間的部分序列在每條read中都存在，可能是潛在的接頭序列。

4.2 clustalw多重序列比對

使用在線clustalw多重序列比對軟件來檢查下我們的top10序列。一般情況下，如果實驗做的沒問題，個別miRNA的表達量會很高，它們可能占據(jù)了絕大部分的reads，因此我們可以利用這個大概率來進行校驗。

圖5. Clustalw輸入top10序列

圖6. 多重序列比對結(jié)果

其中星號為所有序列中都一樣的堿基?？梢钥闯?，后邊的星號連續(xù)排列，因此初步判斷為接頭序列，而星號前面一個堿基僅在test_8里邊與其他的9條不一樣，因此，我們也把T作為接頭的一部分，于是預(yù)測的接頭序列就是TGGAATTCTCGGGTGCCAAGGAACTCC……，經(jīng)過比較，發(fā)現(xiàn)與我們上面提到的接頭序列一樣，因此，可以完全確定接頭序列是正確的。

5，blat比對驗證read組成

在找到接頭序列后，我們可以將接頭序列及以后的堿基全去掉，這樣我們就得到了可能包含miRNA的reads。

圖7. 包含接頭的reads

如圖7所示，綠色為接頭序列，其左側(cè)為可能的miRNA序列，將這10條去接頭的序列（圖8），粘貼到UCSC genome browser的blat工具中，檢索這些序列。

圖8. 去掉接頭的序列

圖9. Blat提交頁面

Blat結(jié)果如下：

圖10. Blat結(jié)果

從圖10的blat結(jié)果可以看出，我們查詢的序列是26或者27 bp，其中identity為100%的序列長度跨了22或者23 bp，因此可以推測序列前或者后加了4 bp的短序列。

打開details鏈接（圖11）查看詳細比對結(jié)果，可以看出都是后邊4 bp的短序列沒有比對上。

圖11. blat詳細比對情況

6，結(jié)論

因此，綜合上面clustalw鑒定的接頭序列，和blat鑒定的read組成，我們得出接過來：該fastq的read組成是圖12的形式。Read結(jié)構(gòu)分析清楚后，我們就可以按部就班地進行后續(xù)的表達、差異、靶基因等分析了。

圖12. Read構(gòu)成

總結(jié)：認清問題，尋找解決問題的思路，思路清晰后，就可以找軟件工具進行處理。

微生信助力發(fā)文章，谷歌引用630+，知網(wǎng)引用470+