RNAseq的第二步就是Aligment, 一般都是有參考基因組的比對竖配，據(jù)同事介紹通今，比對這一步，比對率再80%以上才算是正常的娇妓，很多樣本的比對結(jié)果都能達到90%以上像鸡。
然而，然而哈恰，然而只估，對于第一次跑流程的我，第一個項目着绷，居然不到50% ；赘啤！蓬戚！
QC 統(tǒng)計蘸炸， 均在90%以上堰怨，說明測序數(shù)據(jù)質(zhì)量尚可
遇到問題：Aligment 統(tǒng)計瞬内，比對率 在 28%-47%之間麻裳。
解決思路：找到比對不上的序列幕垦，進行blast 比對循诉，查看比對上的是什么渡嚣。
具體排查：
Aligment 下會有 比對上的fq ID; (A)
QC 中有 所有fq ID （B）
QC 中不含Aligment_fq_ID 的剩下所有ID 為 比對不上的ID （B-A）根據(jù)相應(yīng)軟件(blast)洁仗，查看這些比對不上的ID 到底是什么原因?qū)е碌摹?/p>
簡單介紹用到的軟件
1. samtools
如今序列比對已成為各種生物學(xué)分析中不可缺少的重要環(huán)節(jié)讲衫，通過將未知的基因片段與已知具體信息的基因或基因組進行比較缕棵，并分析其中的相同部分與差異部分孵班，就可以得到該基因片段SNP位點、所屬物種以及可能具有的生物學(xué)功能等重要信息招驴。sam與bam是兩種最常用的比對結(jié)果輸出文件格式篙程，（如轉(zhuǎn)錄組STAR分析軟件輸出的比對結(jié)果為.bam文件等）
bam文件格式是sam文件的二進制格式，占用的存儲空間更小别厘，更利于節(jié)省存儲資源虱饿，而且bam文件的計算處理也更快，但二進制無法直接查看触趴，這就需要借助于工具查看了氮发。
- Samtools 軟件就是用于處理sam與bam格式的工具軟件，能夠?qū)崿F(xiàn)二進制查看冗懦、格式轉(zhuǎn)換爽冕、排序及合并等功能，結(jié)合sam格式中的flag披蕉、tag等信息颈畸，還可以完成比對結(jié)果的統(tǒng)計匯總。同時利用linux中的grep嚣艇、awk等操作命令承冰，還可以大大擴展samtools的使用范圍與功能
基本命令
samtools view [options] <輸入bam文件>
- 主要參數(shù):
-b：以bam格式輸出
-u：以未壓縮的bam格式輸出，一般與linux命令配合使用
-h：在sam輸出中包含header信息
-H：只輸出header信息
-f [INT]：只輸出在比對flag中包含該整數(shù)的序列信息
-F [INT]：跳過比對flag中含有該整數(shù)的序列
-o [file]：標準輸出結(jié)果文件
bam/sam文件 每一列的內(nèi)容
- read名稱食零，通常包括測序平臺等信息
- SAM標記（Flag）困乒，沒有mapping的標記為“ * ”
- chromosome
- 比對上的位置，注意是從1開始計數(shù)贰谣。
- MAPQ（mapping quality娜搂，描述比對的質(zhì)量，數(shù)字越大吱抚，特異性越高百宇，說明該read比對到參考基因組上的位置越唯一）
- CIGAR字串，記錄插入秘豹，刪除携御，錯配以及splice junctions（后剪切拼接的接頭）
- mate名稱，記錄mate pair信息
- mate的位置
- 模板的長度
- read序列
- read質(zhì)量
- 程序用標記
fastq文件格式
FASTQ文件中每個序列通常有四行：
- 序列標識以及相關(guān)的描述信息既绕，以‘@’開頭啄刹；
- 第二行是序列
- 第三行以‘+’開頭，后面是序列標示符凄贩、描述信息誓军，或者什么也不加
- 第四行，是質(zhì)量信息疲扎，和第二行的序列相對應(yīng)昵时，每一個序列都有一個質(zhì)量評分捷雕，根據(jù)評分體系的不同，每個字符的含義表示的數(shù)字也不相同壹甥。
fasta文件格式
在生物信息學(xué)中救巷，F(xiàn)ASTA格式是一種基于文本的、用于表示核苷酸序列或氨基酸序列的格式盹廷。
FASTA文件以 序列表示 和 序列 作為一個基本單元征绸，各行記錄信息如下：
第一行是由大于號">"開頭的任意文字說明。用于序列標記俄占，為了保證后續(xù)分析軟件能夠區(qū)分每條序列管怠，單個序列的標識必須具有唯一性；
-
第二行開始為序列本身缸榄。只允許使用既定的核苷酸或氨基酸編碼符號渤弛。通常核苷酸符號大小寫均可，而氨基酸常用大寫字母甚带。使用時應(yīng)注意有些程序?qū)Υ笮懹忻鞔_要求她肯。文件每行的字母一般不應(yīng)超過80個字符。
fasta 將FASTA序列與質(zhì)量數(shù)據(jù)放到一起鹰贵，就變成了FASTQ
去掉FASTQ文件的后兩行=FASTA
2. blast
- blastn: 核酸序列比對
- blastp: 蛋白質(zhì)序列比對
- blastx: 核酸翻譯成蛋白質(zhì)序列晴氨，再與蛋白質(zhì)序列庫比對
- tblastn:核酸序列庫翻譯成蛋白質(zhì)序列庫，再與蛋白質(zhì)序列比對
- tblastx:核酸與核酸序列庫都翻譯為蛋白質(zhì)序列碉输，再在蛋白質(zhì)水平上比對
- makeblastdb:建立自定義的比對序列庫
使用方法
安裝路徑\bin\blastn -query 比對序列文件名.fasta -db database -out 比對結(jié)果文件名 -evalue le-5 -outfmt 0(或其他比對結(jié)果參數(shù))
- 可以看到籽前，blastn 要求的是fasta格式的數(shù)據(jù)。
實現(xiàn)具體步驟 (eg. pLVX-ShRNA-1)
- 提取Aligment之后（bam文件）匹配上的fq ID
先讀取Align.sorted.bam 文件,提取bam第一列（f1）fq ID列敷钾，并保存到Paired_ID_list.txt文件里
- samtool view + bam文件
$ samtools view Alignment/pLVX-ShRNA-1/Align.sorted.bam | cut -f1 >pLVX-ShRNA-1_unpaired_blast/Paired_ID_list.txt
Paired_ID_list
>>> import pandas as pd
>>> import numpy as np
>>> id_list = pd.read_table("Paired_ID_list.txt",sep = "\t",header=None)
>>> id_list.head()
0
0 A00202:784:HKJKKDSXY:3:1138:3369:12305
1 A00202:784:HKJKKDSXY:3:1123:15953:30968
2 A00202:784:HKJKKDSXY:3:2241:6894:18818
3 A00202:784:HKJKKDSXY:3:1138:3369:12305
4 A00202:784:HKJKKDSXY:3:2241:6894:18818
>>> unilist = np.unique(id_list.values)
>>> print("matchID unique number = ", len(unilist))
matchID unique number = 11704064
-
提取QC之后的所有的fq ID
trim.out
before QC
after QC
$ less ./QC/pLVX-ShRNA-1/pLVX-ShRNA-1_val_1.fq.gz
$ less ./QC/pLVX-ShRNA-1/pLVX-ShRNA-1_val_2.fq.gz
為了測試枝哄，這里只提取前1000條fastq數(shù)據(jù)，
$ zcat QC/pLVX-ShRNA-1/pLVX-ShRNA-1_val_1.fq.gz | head -n 1000 > pLVX-ShRNA-1_unpaired_blast/pLVX-ShRNA-1_val_1_tmp
$ zcat QC/pLVX-ShRNA-1/pLVX-ShRNA-1_val_2.fq.gz | head -n 1000 > pLVX-ShRNA-1_unpaired_blast/pLVX-ShRNA-1_val_2_tmp
檢查一下行數(shù)
$ wc -l pLVX-ShRNA-1_val_1_tmp
1000 pLVX-ShRNA-1_val_1_tmp
下面就是進行ID匹配了阻荒。
pLVX-ShRNA-1_val_1_tmp中仍舊是fastq的標準4行的格式挠锥，我們需要提取其中的seqID，然后檢查該seqID是否在上面的Paired_ID_list.txt當中侨赡，我們需要提取不在的部分fastq序列進行分析蓖租。
- 提取fastq的seqID
方法1：根據(jù)fastq標準格式的規(guī)律，提取seqID
fastqidlst = []
with open("pLVX-ShRNA-1_val_1_tmp",'r') as fastq:
for line in fastq:
if line.startswith('@'):
fastqid = line.strip().split()[0][1:]
fastqidlst.append(fastqid)
方法2：借助于 python Bio::SeqIO
from Bio import SeqIO
fastqidlst = []
for record in SeqIO.parse("pLVX-ShRNA-1_val_1_tmp",format="fastq"):
fastqidlst.append(record.id)
但是羊壹，我們需要的是unpaired seqID對應(yīng)的fastq序列蓖宦。
seqIO可以一步到位，篩選seqID, 并保存seqID對應(yīng)的fastq文件
## 根據(jù) id 來獲炔耙础（過濾）fastq
fq1 = open("pLVX-ShRNA-1_val_1_tmp", "r")
fq2 = open("pLVX-ShRNA-1_val_2_tmp", "r")
OUT1=open("pLVX-ShRNA-1_val_1_tmp_unpaired1","w")
OUT2=open("pLVX-ShRNA-1_val_2_tmp_unpaired2","w")
for record in SeqIO.parse(fq1,format="fastq"):
if record.id not in unilist:
OUT1.write(record.format('fastq'))
for record in SeqIO.parse(fq2,format="fastq"):
if record.id not in unilist:
OUT2.write(record.format('fastq'))
查看fastq文件
less -S pLVX-ShRNA-1_val_2_tmp_unpaired2
pLVX-ShRNA-1_val_2_tmp有 1000條數(shù)據(jù)球昨，居然有496條沒匹配上尔店。
wc -l pLVX-ShRNA-1_val_2_tmp_unpaired2
496 pLVX-ShRNA-1_val_2_tmp_unpaired2
- 利用blastn比對unpaired fastq
blastn輸入的是fasta格式眨攘，首先需要將fastq轉(zhuǎn)成fasta
awk '{if(NR%4 == 1){print ">" substr($0, 2)}}{if(NR%4 == 2){print}}' pLVX-ShRNA-1_val_2_tmp_unpaired2 > pLVX-ShRNA-1_val_2_tmp_unpaired2.fasta
檢查一下主慰，看是不是fasta格式的文件
$ less -S pLVX-ShRNA-1_val_2_tmp_unpaired2.fasta
>A00202:784:HKJKKDSXY:3:1101:5412:1000 2:N:0:GACTAGTA+TCTTTCCC
TGTCGCTCCATCAAGCTTTCGCTCATTGTGGAAAATTCCTTACTGCTGCCTCCCGTAGGAGTCTGGGCCGTATCTCAGTCCCAGTGTGGCCGTACAGCCTCTCGGCCCGGCTAAACATCATCGTCTTGGTAGGCCATTACCCTAC
blastn進行核酸序列比對
/home/Software/ncbi-blast-2.11.0+/bin/blastn -task blastn -query pLVX-ShRNA-1_val_2_tmp_unpaired2.fasta \
-db /home/database/nt_blast/nt -show_gis -evalue 1e-3 -max_target_seqs 5 -num_threads 6 -out unpaired2_blast.nt.aln \
-outfmt "6 qseqid qlen sseqid slen qstart qend sstart send length pident nident bitscore evalue staxids stitle"
發(fā)現(xiàn)都是支原體污染，問了一下養(yǎng)過細胞的小伙伴鲫售，說是細胞被支原體污染是件很常見的事情啊啊啊啊啊啊啊啊啊啊啊啊啊
