前段时间送去测序的实验数据拿回来了，然后做完比对之后发现unique_reads比对率不够理想，确保分析数据过程中没有出错的前提下，个人认为实验过程中可能有污染。所以需要随机抽取未比对上的reads进行BLAST寻找出可能的污染源。

一、重温bam文件格式是关键！

        bam文件有12列：举个例子才能更直观。(哈哈哈，其实是我自己也记不住每列代表的啥意思，所以一般都会先查看一下bam文件然后再直接看每列代表的啥意思）

        整整齐齐的查看bam文件：samtools view x.bam | less -S

第一列：QNAME:比对序列的名称。

第二列：FLAG:比对的类型。paring、strand、 mate strand等等。不同的数值代表不一样的比对类型。

第三列：RNAME：表示read比对的那条序列的序列名称。如果第三列是”*“，则说明这条read没有比对上。

第四列：POS。表示read比对到RNAME这条序列最左边的位置。

第五列：mapping的质量值。

第六列：CIGAR。比对的具体情况。

第七列：MRNM。mate序列所在参考序列的名称。mate一般指大的片段序列。

第八列：MPOS。mate序列在参考序列上的位置。

第九列：ISIZE。文库插入片段长度。

第十列：reads sequence。

第十一列：reads对应的ASCII码格式的碱基质量值。

第十二列：可选的区域。

---

二、从bam文件中提取未比对上的reads。

       接下来要做的就是把bam文件中的第三列是"*"对应的第十列的序列找出来就可以啦，哈哈哈，我真是个小机灵！

      samtools view x.bam | awk '$3=="*" {print ">"$1"\n"$10}' >x_no_mapped_reads.txt

什么是快乐星球？！

目前就是轻松得到了未比对上的reads 序列！

嘻嘻嘻

---

三、NCBI进行BLAST

进入NCBI官网点击BLAST 点击Nucleotide BLAST(如果你是其他序列就点击对应的BLAST工具） 随机挑选某条未比对上的reads序列，然后点击show results in new window 最后点击BLAST就可以啦 大功告成！