? ? ? ?前段時間送去測序的實驗數(shù)據(jù)拿回來了榜轿，然后做完比對之后發(fā)現(xiàn)unique_reads比對率不夠理想汉规，確保分析數(shù)據(jù)過程中沒有出錯的前提下包斑，個人認(rèn)為實驗過程中可能有污染摘昌。所以需要隨機(jī)抽取未比對上的reads進(jìn)行BLAST尋找出可能的污染源斥杜。
一虱颗、重溫bam文件格式是關(guān)鍵沥匈！
? ? ? ? bam文件有12列：舉個例子才能更直觀。(哈哈哈忘渔，其實是我自己也記不住每列代表的啥意思高帖，所以一般都會先查看一下bam文件然后再直接看每列代表的啥意思）
? ? ? ? 整整齊齊的查看bam文件：samtools view x.bam | less -S
第一列：QNAME:比對序列的名稱。
第二列：FLAG:比對的類型畦粮。paring散址、strand、 mate strand等等宣赔。不同的數(shù)值代表不一樣的比對類型预麸。
第三列：RNAME：表示read比對的那條序列的序列名稱。如果第三列是”*“儒将，則說明這條read沒有比對上吏祸。
第四列：POS。表示read比對到RNAME這條序列最左邊的位置椅棺。
第五列：mapping的質(zhì)量值犁罩。
第六列：CIGAR。比對的具體情況两疚。
第七列：MRNM床估。mate序列所在參考序列的名稱。mate一般指大的片段序列诱渤。
第八列：MPOS丐巫。mate序列在參考序列上的位置。
第九列：ISIZE勺美。文庫插入片段長度递胧。
第十列：reads sequence。
第十一列：reads對應(yīng)的ASCII碼格式的堿基質(zhì)量值赡茸。
第十二列：可選的區(qū)域缎脾。
二、從bam文件中提取未比對上的reads占卧。
? ? ? ?接下來要做的就是把bam文件中的第三列是"*"對應(yīng)的第十列的序列找出來就可以啦遗菠，哈哈哈，我真是個小機(jī)靈华蜒！
? ? ? samtools view x.bam | awk '$3=="*" {print ">"$1"\n"$10}' >x_no_mapped_reads.txt
什么是快樂星球辙纬？！
目前就是輕松得到了未比對上的reads 序列叭喜！
嘻嘻嘻
三贺拣、NCBI進(jìn)行BLAST
進(jìn)入NCBI官網(wǎng)點擊BLAST
點擊Nucleotide BLAST(如果你是其他序列就點擊對應(yīng)的BLAST工具）
隨機(jī)挑選某條未比對上的reads序列，然后點擊show results in new window 最后點擊BLAST就可以啦
大功告成！
