参考：
生信：2：sam格式文件解读
FastQ Screen -- 调查fastq测序数据是否污染
比对率不理想-污染检测
这或许是我写的最全的BLAST教程
linux下安装blast并创建nt数据库
什么！！！超70G的NT数据库文件一个小时搞定？
假如Fastqc中GC% 出现双峰怎么判断是否可能是其他物种混入呢？
我们可以提取序列进行blast
image.png
1.blast 本地比较
- 安装seqtk
git clone https://github.com/lh3/seqtk.git;
cd seqtk; make
- 安装blast+ 参考：https://zhuanlan.zhihu.com/p/141644734
wget https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.11.0+-x64-linux.tar.gz
tar -zxvf ncbi-blast-2.11.0+-x64-linux.tar.gz
vi ~/.bashrc
export PATH=/youpath/ncbi-blast-2.11.0+/bin:$PATH
- 下载nt/nr库（核酸/蛋白质）
也可以考虑自己构建库
(base) [11:28:39] kcao@login:~/genome_human/ncbi
$ ascp -v -k 1 -T -l 200m -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh anonftp@ftp.ncbi.nlm.nih.gov:/blast/db/FASTA/nt.gz ./
$ ascp -v -k 1 -T -l 200m -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh anonftp@ftp.ncbi.nlm.nih.gov:/blast/db/FASTA/nr.gz ./
- 运行blastn 查看哪些基因组比对较多
### 1.提取测试fa
~/tools/seqtk/seqtk sample -s100 ACT-08_trimmed.R1.fq.gz 100000 |~/tools/seqtk/seqtk seq -a - >Sample.20W.fa &&
~/tools/seqtk/seqtk sample -s100 ACT-08_trimmed.R2.fq.gz 100000 |~/tools/seqtk/seqtk seq -a - >>Sample.20W.fa
### 2.blastn 比对
$ module load BLAST+/2.6.0-foss-2016b-Python-2.7.11
blastn -query Sample.20W.fa -db 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
less Sample.20W.nt.txt|awk '{a[$1]++}END{for( i in a){print i,a[i] | "sort -nrk 2"}}' |head
2.网页版blast
假如比对率只有50% 左右，同时Fastqc 出现了双峰，可能怀疑有污染。
2.1 提取未比对的行
head -100000 input-test.sam|awk '($3=="*"){print $0}'|head
image.png
2.2 提取前20行没有比对的reads序列。
cat input-EAF1.sam|awk '($3=="*"){print $10}'|head -20
image.png
2.3 转换成fa格式
cat input-test.sam|awk '($3=="*"){n++;print ">unmapping"n"\n"$10}'|head -40
image.png
2.4 将fasta文件上传到blast中比对，找到可疑的污染。(可能玉米样本混入)
image.jpeg
2.5 比对到玉米基因组上
[kcao@login sam_file]$
echo "bowtie2 -p 8 -x /public/home/kcao/Zea_mays.AGPv3.30/Zea_mays.AGPv3.30 -1 /public/home/kcao/tools/CSATK/output/04_trim/input-test_1.fastq -2 /public/home/kcao/tools/CSATK/output/04_trim/input-test_2.fastq -S /public/home/kcao/tools/CSATK/output/bowtie2_mapping/input-EAF1.mapping.Zea.may.sam"|qsub -d ./ -N Zea.may_input_test.mapping
echo "samtools view -bS -o input-test.mapping.Zea.may.sam.bam input-test.mapping.Zea.may.sam "|qsub -d ./ -N samtobam
samtools flagstat input-test.mapping.Zea.may.sam.bam >input-test.mapping.Zea.may.sam.bam.flagstat"|qsub -d ./ -N flagstat
统计结果表明：发现比对到玉米上比较多
image.png
3.FastQ Screen 工具
好像也要指定可能的污染源
image.png
思考
- 本地blast可能效果更好，可以查看比对到不同物种比例，但是nt/nr数据库非常大，建库或者下载耗时很长。
- 也可以考虑将常见的物种索引建好，建一个在线服务website，提取部分序列进行BWA，看比对率。
欢迎评论交流~
网友评论