1.从bam文件中提取未比对上的reads到新bam文件
samtools view -bf 4 IP_988_1_H3K4me3.bam > IP_988_1_H3K4me3unmapped.bam
2.提取出未比对上的双端reads （bamtofastq）
bamToFastq -i IP_988_1_H3K4me3.bam -fq IP_988_1_H3K4me3unmappedR1.fastq -fq2 IP_988_1_H3K4me3unmappedR2.fastq
3.从双端各随机抽取1W序列，组成2W条query序列（seqtk seq -a为转化fastq到fa）。
seqtk sample -s100 Sample.R1.fastq.gz 10000 |seqtk seq -a - >Sample.2W.fa &&
seqtk sample -s100 Sample.R2.fastq.gz 10000 |seqtk seq -a - >>Sample.2W.fa
4.下载NT和NR数据库
ascp -v -k 1 -T -l 200m -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh （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 ./
注：需要提前下载aspera-cli；该软件可以从大安装；安装命令为：conda install -y -c hcc aspera-cli
5.数据库比对
gunzip nt.gz
makeblastdb -in nt -out nt -dbtype nucl
blastn -query Sample.2W.fa -db nt -num_threads 8 -outfmt '6 staxids qseqid sseqid pident length mismathch gapopen qstart qend sstart send evalue bitscore qcovs' -evalue 1e-10 -max_target_seqs 1 -out Sample.2W.nt.txt
6.提取比对结果
利用awk整理比对结果，并统计比对到各物种的比例
初步确定存在病菌类和人源污染，去找公司argume！
