gmap使用方法

【提示】本文写作命令为：abysw blog gmap使用方法

有很多软件需要反复使用，但是总是会忘了怎么用，需要备注一下。

Consensus sequences, FLNC sequences, and nonredundant transcripts were mapped to the Z. japonica reference transcript sequences (http:// zoysia. kazusa. or. jp) using GMAP [50].
Then, PacBio sequencing was mapped to the bamboo genome using GMAP with the following option: -B 5 -K 8000 -t 40 -f 2 -n 1 (Wu et al., 2016).

安装gmap

在maker安装的时候已经安装好了，可以直接用

1. 建库
   gmap_build -D ./ -d ${S}_genome ${S}_genome.fa

2. 比对
   gmap -B 5 -K 8000 -t 40 -f 2 -n 1 -D ./ -d ${S}_genome $R.transcripts.hq.fasta > $R.transcripts.hq.fasta.gmap 2> $R.transcripts.hq.fasta.fail

这里我测试的比对率是98.85%。

看看参数：
-B, –batch=INT Batch mode (default = 2)
-t, –nthreads=INT Number of worker threads
-f, –format=INT Other format for output (also note the -A and -S options
-n, –npaths=INT Maximum number of paths to show (default 5). If set to 1, GMAP
will not report chimeric alignments, since those imply
two paths. If you want a single alignment plus chimeric
alignments, then set this to be 0.
上面的三个参数没什么用，-n参数可以调试一下，看看效果。
–min-trimmed-coverage=FLOAT Do not print alignments with trimmed coverage less this value (default=0.0, which means no filtering)
–min-identity=FLOAT Do not print alignments with identity less this value (default=0.0, which means no filtering)
也可以调试这两个参数的用法。

-K 
--max-intronlength-middle=INT  Max length for one internal intron (default 500000).  Note: for backward
                               compatibility, the -K or --intronlength flag will set both
                               --max-intronlength-middle and --max-intronlength-ends.
                               Also see --split-large-introns below.
--max-intronlength-ends=INT    Max length for first or last intron (default 10000).  Note: for backward
                               compatibility, the -K or --intronlength flag will set both
                               --max-intronlength-middle and --max-intronlength-ends.

================================== 以下为调试区 ==============================
gmap -B 5 -K 8000 -t 40 -f 2 -n 0 -D ./ -d ${S}_genome $R.transcripts.hq.fasta > $R.transcripts.hq.fasta.gmap0ic 2> $R.transcripts.hq.fasta.gmap.fail0ic &
n=0 多输出202个嵌合体。

gmap -B 5 -K 8000 -t 40 -f 2 -n 0 –min-trimmed-coverage=0.85 –min-identity=0.9 -D ./ -d ${S}_genome $R.transcripts.hq.fasta > $R.transcripts.hq.fasta.gmap0ic 2> $R.transcripts.hq.fasta.gmap.fail0ic &
覆盖度和相似度过滤的结果：
未匹配序列：856 –> 1853， 比对率为97.52。
嵌合体序列：170条

gmap –min-trimmed-coverage=0.85 –min-identity=0.9 -D ./ -d ${S}_genome $R.transcripts.hq.fasta > $R.transcripts.hq.fasta.gmapic 2> $R.transcripts.hq.fasta.gmap.failic &

测试最基本的参数：
gmap -B 5 -t 40 -f 2 -n 1 -D ./ -d ${S}_genome $R.transcripts.hq.fasta > $R.transcripts.hq.fasta1.gff 2> $R.transcripts.hq.fasta1.gmap.err &
gmap -B 5 -t 40 -f 2 -n 0 -D ./ -d ${S}_genome $R.transcripts.hq.fasta > $R.transcripts.hq.fasta.gff 2> $R.transcripts.hq.fasta.gmap.err &
两者的比对率都是98.85%；
嵌合体为213个（较之前多11个嵌合体，为-K参数所致）。
path2 暂时不知道有没有用，先用文献中的参数吧。
部分基因跨度太大。
==============================================================================
gmap_build -D ./ -d ${S}_genome ${S}_genome.fa
gmap -B 5 -K 8000 -t 40 -f 2 -n 1 -D ./ -d ${S}_genome $R.transcripts.hq.fasta > $R.transcripts.hq.fasta.gff 2> $R.transcripts.hq.fasta.gmap.err &

#################文件默认变量#####################
$G genome fasta file
$S species
$R NGS reads / NGS sample name
##################################################