一湿痢、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序列比對(duì)到參考序列上购桑，
${sample}_contig.fa: target 參考序列群嗤，可以是Megahit組裝后的contigs或者參考基因組${sample}_R[12].fq: query 查詢序列菠隆，paired-end則要指定R1和R2-t： Number of threads [3].-a： 生成CIGAR, 并以SAM格式輸出比對(duì)結(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)比對(duì)到參考基因組 -
-x sr： Short single-end reads without splicing.
-
-o FILEOutput alignments to FILE [stdout].
$minimap2 -ax sr \ # 短reads比對(duì)模式；以SAM格式輸出比對(duì)結(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比對(duì)結(jié)果文件
二、samtools view
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 FILEoutput file name [stdout]
1. 從bam中提取mapped/unmapped 的序列信息
- 從比對(duì)結(jié)果文件
${sample}_aln.bam中提取匹配上的序列信息誉碴，并保存到${sample}_mapped.bam文件
${samtools} view -bF 4 ${sample}_aln.bam > ${sample}_mapped.bam
- 從比對(duì)結(jié)果文件
${sample}_aln.bam中提取沒有匹配上的序列信息，并保存到${sample}_unmapped.bam文件
${samtools} view -bf 4 ${sample}_aln.bam > ${sample}_unmapped.bam
- 從比對(duì)結(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是序列比對(duì)厚度輸出格式管搪，包括頭部信息（以@開頭）和比對(duì)信息兩部分組成
- Header 信息
-
@HD VN:1.0 SO:unsorted【排序類型】
頭部區(qū)第一行：VN是格式版本更鲁；SO表示比對(duì)排序的類型澡为，有unknown（default）媒至，unsorted驯绎，queryname和coordinate幾種谋旦。samtools軟件在進(jìn)行行排序后不能自動(dòng)更新bam文件的SO值册着，而picard卻可以乞巧。 -
@SQ SN:contig1 LN:9401【參考序列ID及其長度】
參考序列名绽媒，這些參考序列決定了比對(duì)結(jié)果sort的順序是辕，SN是參考序列名；LN是參考序列長度疙教；每個(gè)參考序列為一行。
例如：@SQ SN:NC_000067.6 LN:195471971 -
@RG ID:sample01【樣品基本信息】
Read Group裸弦。1個(gè)sample的測(cè)序結(jié)果為1個(gè)Read Group理疙；該sample可以有多個(gè)library的測(cè)序結(jié)果，可以利用bwa mem -R 加上去這些信息。
例如：@RG ID:ZX1_ID SM:ZX1 LB:PE400 PU:Illumina PL:Miseq
ID：樣品的ID號(hào) SM：樣品名 LB：文庫名 PU：測(cè)序以 PL：測(cè)序平臺(tái)
這些信息可以在形成sam文件時(shí)加入粤攒，ID是必須要有的后面是否添加看分析要求 -
@PG ID:bowtie2 PN:bowtie2 VN:2.0.0-beta7【比對(duì)所使用的軟件及版本】
例如：@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格式文件解讀
3. 從上面minimap2的比對(duì)輸出文件${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中.
其對(duì)序列的分類依據(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.
對(duì)于PE測(cè)序reads, 同時(shí)指定-1 R1.fq -2 R2.fq -s singleton.fq時(shí)，samtools會(huì)將 flag1和flag2配對(duì)的序列分別輸出到-1 -2指定的文件骗露，對(duì)于無法匹配的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 -
