minimap2 + samtools 比對參考序列，并提取unmapped reads

一、minimap2

Manual Reference Pages - minimap2 (1)

  • Long-read alignment with CIGAR:

minimap2 -a [-x preset] target.mmi query.fa > output.sam

minimap2 -c [-H] [-k kmer] [-w miniWinSize] [...] target.fa query.fa > output.paf

minimap2 將query序列比對到參考序列上，

  • ${sample}_contig.fa : target 參考序列，可以是Megahit組裝后的contigs或者參考基因組

  • ${sample}_R[12].fq: query 查詢序列，paired-end則要指定R1和R2

  • -t： Number of threads [3].

  • -a ： 生成CIGAR, 并以SAM格式輸出比對結(jié)果（minimap2默認(rèn)輸出PAF格式的文件）

  • -x [STR]： 預(yù)設(shè)選項(xiàng)。部分選項(xiàng)：

    • -x map-ont： 默認(rèn)選項(xiàng)， 將noisy long reads(~10% error rate)比對到參考基因組
    • -x sr： Short single-end reads without splicing.

  • -o FILE Output alignments to FILE [stdout].

$minimap2 -ax sr \          # 短reads比對模式；以SAM格式輸出比對結(jié)果
    -t ${threads} \         # 設(shè)置CPU線程數(shù)（threads >= 3）
    3_assembly/${sample}_contig.fa \    # target.fa
    ${sample}_R1.fq \          # query.R1.fq
    ${sample}_R2.fq \          # query.R2.fq
    -o ${sample}_aln.sam                 # minimap2比對結(jié)果文件

二、samtools view

  1. Manual page from samtools-1.15
  2. SAM格式詳細(xì)說明-簡書
  3. samtools常用命令詳解

sam轉(zhuǎn)bam格式

# 將sam文件轉(zhuǎn)換成bam文件
${samtools} view -bS abc.sam > abc.bam
# 等價(jià)于
${samtools} view -b -S abc.sam -o abc.bam
  • -b： output BAM。默認(rèn)下輸出是 SAM 格式文件，該參數(shù)設(shè)置輸出 BAM 格式
  • -S： input is SAM。默認(rèn)下輸入是 BAM 文件，若是輸入是 SAM 文件，則最好加該參數(shù)，否則有時(shí)候會(huì)報(bào)錯(cuò)。
  • -o FILE output file name [stdout]

1. 從bam中提取mapped/unmapped 的序列信息

  • 從比對結(jié)果文件${sample}_aln.bam中提取匹配上的序列信息，并保存到${sample}_mapped.bam文件
${samtools} view -bF 4 ${sample}_aln.bam > ${sample}_mapped.bam
  • 從比對結(jié)果文件${sample}_aln.bam中提取沒有匹配上的序列信息，并保存到${sample}_unmapped.bam文件
${samtools} view -bf 4 ${sample}_aln.bam > ${sample}_unmapped.bam
  • 從比對結(jié)果文件${sample}_aln.bam中提取沒有 read1和 read2均匹配上的序列信息，并保存到${sample}_all.mapped.bam文件
${samtools} view -bF 12 ${sample}_aln.bam > ${sample}_all.mapped.bam

2. 給予指定的reference文件，將sam轉(zhuǎn)化為bam

sam是序列比對厚度輸出格式，包括頭部信息（以@開頭）和比對信息兩部分組成

  • Header 信息
    • @HD VN:1.0 SO:unsorted 【排序類型】
      頭部區(qū)第一行：VN是格式版本；SO表示比對排序的類型，有unknown（default），unsorted，queryname和coordinate幾種。samtools軟件在進(jìn)行行排序后不能自動(dòng)更新bam文件的SO值，而picard卻可以。
    • @SQ SN:contig1 LN:9401 【參考序列ID及其長度】
      參考序列名，這些參考序列決定了比對結(jié)果sort的順序，SN是參考序列名；LN是參考序列長度；每個(gè)參考序列為一行。
      例如：@SQ SN:NC_000067.6 LN:195471971
    • @RG ID:sample01 【樣品基本信息】
      Read Group。1個(gè)sample的測序結(jié)果為1個(gè)Read Group；該sample可以有多個(gè)library的測序結(jié)果，可以利用bwa mem -R 加上去這些信息。
      例如：@RG ID:ZX1_ID SM:ZX1 LB:PE400 PU:Illumina PL:Miseq
      ID：樣品的ID號(hào) SM：樣品名 LB：文庫名 PU：測序以 PL：測序平臺(tái)
      這些信息可以在形成sam文件時(shí)加入，ID是必須要有的后面是否添加看分析要求
    • @PG ID:bowtie2 PN:bowtie2 VN:2.0.0-beta7 【比對所使用的軟件及版本】
      例如：@PG ID:bwa PN:bwa VN:0.7.12-r1039 CL:bwa sampe -a 400 -f ZX1.sam -r @RG ID:ZX1_ID SM:ZX1 LB:PE400 PU:Illumina PL:Miseq …/0_Reference/Reference_Sequence.fa ZX_HQ_clean_R1.fq.sai ZX_HQ_clean_R2.fq.sai …/2_HQData/ZX_HQ_clean_R1.fq …/2_HQData/ZX_HQ_clean_R2.fq
      這里的ID是bwa，PN是bwa，VN是0.7.12-r1039版本。CL可以認(rèn)為是運(yùn)行程序@RG是上面RG表示的內(nèi)容，后面是程序內(nèi)容，這里的@GR內(nèi)容是可以自己在運(yùn)行程序是加入的

sam轉(zhuǎn)bam

(1) sam文件的header包含@SQ ，即sam中包含了reference的信息。

${samtools} view -bS aln.sam > aln.bam

(2) sam文件不包含header或者h(yuǎn)eader不包含@SQ ，即sam中不包含了reference的信息，此時(shí)需要提供生成sam文件時(shí)使用的reference文件。

${samtools} faidx ref.fa  # 建索引
${samtools} view -bS -t ref.fa.fai aln.sam > aln.bam  # 輸出轉(zhuǎn)換結(jié)果
# 或者
${samtools} view -bS -T ref.fa aln.sam > aln.bam # 自動(dòng)建索引，并輸出轉(zhuǎn)換結(jié)果

ref: 1. 生信：2：sam格式文件解讀

  1. Manual page from samtools-1.15: samtools-view

3. 從上面minimap2的比對輸出文件${sample}_aln.sam文件中，提取未匹配的文件，并保存為bam文件

${samtools} view -bS \
    -T 3_assembly/${sample}_contig.fa \
    -f 4 \
    ${sample}_aln.sam
    > ${sample}_aln.bam

三、samtools fastq

samtools fastq [options] in.bam
將輸入的bam文件轉(zhuǎn)化為fastq文件，并將結(jié)果保存至指定的輸出-1 -2 -o -0-s中.
其對序列的分類依據(jù)是序列末尾的READ1 flag和READ2 flag, flag有3類：

  • 1 : Only READ1 is set.
  • 2 : Only READ2 is set.
  • 0 : Either both READ1 and READ2 are set; or neither is set.

對于PE測序reads, 同時(shí)指定-1 R1.fq -2 R2.fq -s singleton.fq時(shí)，samtools會(huì)將 flag1和flag2配對的序列分別輸出到-1 -2指定的文件，對于無法匹配的flag 1/2，全部的flag 0 都會(huì)保存到-s 指定的文件中。

refs: Manual page from samtools-1.15: samtools fastq

${samtools} fastq -n \
    -1 3_assembly/${sample}_unmapped_R1.fq.gz \
    -2 3_assembly/${sample}_unmapped_R2.fq.gz \
    -s 3_assembly/${sample}_unmapped_singleton.fq.gz \
    ${sample}_aln.bam

四、samtools 給予管道（|）的輸入輸出 -

samtools旨在處理數(shù)據(jù)流（stream），
它可以將-作為管道中的標(biāo)準(zhǔn)輸入文件（stdin）;
也可以將- 作為管道中的標(biāo)準(zhǔn)輸出文件（stdout）.

全部代碼

threads=12
sample=A01
minimap2=/software/miniconda2/envs/metadenovo/bin/minimap2
samtools=/software/samtools/samtools-1.8/bin/samtools
## mapping
$minimap2 -ax sr -t ${threads} \
    3_assembly/${sample}_contig.fa \ 
    ${sample}_R1.fq \
    ${sample}_R2.fq | \
${samtools} view -bS -T 3_assembly/${sample}_contig.fa \
    -f 4 - | \
${samtools} fastq -n \
    -1 3_assembly/${sample}_unmapped_R1.fq.gz \
    -2 3_assembly/${sample}_unmapped_R2.fq.gz \
    -s 3_assembly/${sample}_unmapped_singleton.fq.gz -
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。
禁止轉(zhuǎn)載，如需轉(zhuǎn)載請通過簡信或評論聯(lián)系作者。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容