比对率不理想-污染检测

污染检测，即通过blast，对样本的序列进行nt总库比对，看样本reads中不同物种的占比情况。

在进行基因组比对时，发现部分或者全部样本比对率低于80%，或者是两样本之间比对率相差超过10%，就需要抽取reads（paired）进行污染检测。

首先，我们选择1-2个比对率较高的作为对照，以及比对率低的样本（样本很多的可选部分样本），抽取测序的双端reads。

seqtk sample -s100 Sample.R1.fastq.gz 100000 |seqtk seq -a - >Sample.20W.fa &&
seqtk sample -s100 Sample.R2.fastq.gz 100000 |seqtk seq -a - >>Sample.20W.fa

而后，通过blast比对结果，看比对率低的样本中是掺杂了什么物种的序列。
各抽取10W 条reads，合并20W reads,其中seqtk seq -a 将fq转fa。

blastn -query Sample.20W.fa -db nt/nt -outfmt '6 staxids qseqid sseqid pident length mismathch gapopen qstart qend sstart send evalue bitscore qcovs' -evalue 1e-10 -max_target_seqs 1 -out pollution_test/Sample.20W.nt.txt -num_threads 10
less Sample.20W.nt.txt|awk '{a[$1]++}END{for( i in a){print i,a[i] | "sort -nrk 2"}}' |head

• 如果是其他物种，且占比超过1%，就说明该样本可能有污染；(注意筛选阈值调整，P值，分值以及reads名称去重)

  
  1.1污染检测结果
  

  可以看出。该样本第二位比对情况有异常。

  
  1.2物种信息

• 如果没有明显的外源物种序列或者外源序列占比低于1%，则该样本可能没有污染。

  
  2.未污染样本

进阶：

• 百分比并不是 比对上reads数目/总reads数目；
• 抽取的序列可能是20W条，5W条或者其他，我们需要看一下比对上多少序列；
  以图2为例，抽取5w序列，比对上人，猩猩，狒狒的总共有近2W条，那么其中有3W未比对到人和近源物种，说明该样本可能比较特殊，需要和相关人员进一步沟通。