未比对上的bam reads 处理

查看

bam文件 flags 4代表read未比对上；

samtools view -f 4 *bam|less -S

提取

提取unmapped序列:

samtools view -b -h -f 4 *.bam > unmapped.bam

-h 文件包含header line;-f,提取;-b,输出为bam格式

###-F参数是过滤的意思（filter）；类似grep -v

samtools view -b -h -F 4 *.bam > mapped.bam

F意思为过滤掉以4为标签的序列

提取双端序列都mapping到参考序列（4+8）的结果

samtools view -bF 12 *.bam > mapped.bam

3.bam2fastq

bamToFastq -bam file_unmapped.bam -fq1 unmappedR1.fastq -fq2 unmappedR2.fastq

---

ref:https://www.jianshu.com/p/471a4fbcdd48

wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nt.gz
wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz
wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nt.gz.md5
wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz.md5
md5sum -c nt.gz.md5
md5sum -c nr.gz.md5
tar -xzvf nr.gz
tar -xzvf nt.gz

mkdir nr_db
mkdir nt_db
makeblastdb -in nr -dbtype prot -title make_nr -parse_seqids -out ./nr_db/nr -logfile make_nr.log
makeblastdb -in nt -dbtype nucl -title make_nt -parse_seqids -out ./nt_db/nt -logfile make_nt.log

或者使用NCBI处理好的db。

wget -c ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt*
wget -c ftp://ftp.ncbi.nlm.nih.gov/blast/db/nr*

blastp：蛋白序列与蛋白库作比对，直接比对蛋白序列的同源性。
blastx：核酸序列与蛋白库作比对，将核酸序列先翻译成蛋白序列，再将其与蛋白库作比对。
-blastn：核酸序列与核酸库的比对，直接比对核酸序列的同源性。
tblastn：蛋白序列对核算库的比对，现将核酸库翻译成蛋白库，再将蛋白序列与翻译后的蛋白库进行比对。
tblastx：核酸与核酸数据库在蛋白质水平比较

awk '{if(NR%4 == 1){print ">" substr($0, 2)}}{if(NR%4 == 2){print}}' test.fastq > test.fasta
blastn -query test.fasta -out test.result -db ./nt_db/nt 
## 输出一条最优比对结果
blastn -query test1.fa -out test1.align -db ./nt_db_NCBI/nt -outfmt 6 -subject_besthit -num_threads 4