STAR 2-pass, picard 到gatk的使用
!!!对 RNA-seq 产出的数据进行变异检测分析，与常规重测序的主要区别就在序列比对这一步，因为 RNA-seq 的数据是来自转录本的，比对到参考基因组需要跨越转录剪切位点，所以 RNA-seq 进行变异检测的重点就在于跨剪切位点的精确序列比对。
STAR相比较其他两款软件有较高的唯一比对率；STAR会将没有paired mapping上的reads都剔除，避免single reads比对到基因组上；并且STAR对lower-quality（包括more soft-clipped和错配碱基）比对有较高的容忍度
构建基因组索引
STAR --runThreadN 25 --runMode genomeGenerate --genomeDir index_dir --genomeFastaFiles genome.fasta --sjdbGTFfile genome.gtf --sjdbOverhang 149
--runThreadN：线程数。
--runMode genomeGenerate：构建基因组索引。
--genomeDir：索引目录。（index_dir一定要是存在的文件夹，需提前建好）
--genomeFastaFiles：基因组文件。
--sjdbGTFfile：基因组注释文件。
--sjdbOverhang：reads长度减1。
实例
STAR --runThreadN 25 --runMode genomeGenerate --genomeDir STAR_index --genomeFastaFiles mm10.fa --sjdbGTFfile mm10_genes.gtf --sjdbOverhang 149
####
STAR --runThreadN 25 --runMode genomeGenerate --genomeDir STAR_index --genomeFastaFiles GRCm39.fa --sjdbGTFfile GRCm39.gtf --sjdbOverhang 149
STAR --runThreadN 25 --runMode genomeGenerate --genomeDir STAR_index_mm39 --genomeFastaFiles mm39.fa --sjdbGTFfile mm39.gtf --sjdbOverhang 149
#若出现报错如下：
[u20111230014@cpu23 GCF_GRCm38.p6]$ more Out.star_index
STAR --runThreadN 25 --runMode genomeGenerate --genomeDir STAR_index --genomeFastaFiles GRCm38.fa --sjdbGTFfile GRCm38.gtf --sjdbOverhang 149
STAR version: 2.7.9a compiled: 2021-05-04T09:43:56-0400 vega:/home/dobin/data/STAR/STARcode/STAR.master/source
Mar 29 22:19:17 ..... started STAR run
Mar 29 22:19:17 ... starting to generate Genome files
EXITING because of FATAL PARAMETER ERROR: limitGenomeGenerateRAM=31000000000is too small for your genome
SOLUTION: please specify --limitGenomeGenerateRAM not less than 33524373088 and make that much RAM available
##解决
一般平台使用slurm时，cpu核数与内存之间要相对平衡，建议1:5，即#SBATCH -c 2; #SBATCH --mem=10G。
将脚本的cpu改大一点，并配上 --limitGenomeGenerateRAM 报错的内存数字
#!/bin/bash
#SBATCH -J Job.star_index
#SBATCH -p dna
#SBATCH -N 1
#SBATCH --mem=80G
#SBATCH --cpus-per-task=40
#SBATCH -t 1-10:00:00
#SBATCH -o Out.star_index
#SBATCH --mail-user=394710725@qq.com
STAR --runThreadN 40 --runMode genomeGenerate --genomeDir STAR_index --genomeFastaFiles GRCm38.fa --sjdbGTFfile GRCm38.gtf --sjdbOverhang 149 --limitGenomeGenerateRAM 335
24373088
STAR 2-pass mode
为了发现更加灵敏的new junction，STAR建议使用2-pass mode，其能增加检测到的new junction数目，使得更多的splices reads能mapping到new junction。因此STAR先用一般参数做一遍mapping，收集检测到的junction信息，然后利用这已经annotated junction来做第二次mapping
双端比对
示例
STAR --twopassMode Basic --quantMode TranscriptomeSAM GeneCounts --runThreadN 6 --genomeDir index_dir --alignIntronMin 20 --alignIntronMax 50000 --outSAMtype BAM SortedByCoordinate --sjdbOverhang 149 --outSAMattrRGline ID:sample SM:sample PL:ILLUMINA --outFilterMismatchNmax 2 --outSJfilterReads Unique --outSAMmultNmax 1 --outFileNamePrefix out_prefix --outSAMmapqUnique 60 --readFilesCommand gunzip -c --readFilesIn seq1.fq.gz seq2.fq.gz
##单端比对
示例
STAR \
--runThreadN 20 \
--genomeDir hg19_STAR_db \
--readFilesIn reads.fq \
--sjdbGTFfile hg19.gtf \
--sjdbOverhang 100 \
--outFileNamePrefix sampleA \
--outSAMtype BAM SortedByCoordinate
##命令参数也很好理解：
--runThreadN ：设置线程数
--genomeDir ：index输出的路径
--genomeFastaFiles ：参考基因组序列
--sjdbGTFfile ：参考基因组注释文件
--sjdbOverhang ：这个是reads长度的最大值减1，默认是100
如果是per-sample(单样本)的2-pass mapping，可以直接用--twopassMode Basic参数将第两步mapping中的make index省去，直接再mapping
操作路径：/home/u20111230014/workspace/genome/STAR_mm10
实例
STAR --runThreadN 25 --genomeDir ./STAR_index --readFilesIn DPP-0_clean.fq --sjdbGTFfile mm10_genes.gtf --sjdbOverhang 100 --outFileNamePrefix DPP-0 --outSAMtype BAM SortedByCoordinate
操作路径：/home/u20111230014/workspace/genome/GRCm39/STAR_GRCm39
实例
STAR --runThreadN 25 --genomeDir ./STAR_index --twopassMode Basic --readFilesIn DPP-0_clean.fq --sjdbGTFfile GRCm39.gtf --sjdbOverhang 149 --outFileNamePrefix DPP-0-GRCm39 --outSAMtype BAM SortedByCoordinate
操作路径： /home/u20111230014/workspace/genome/mm10_to_mm39
实例
STAR --runThreadN 25 --genomeDir ./STAR_index --twopassMode Basic --readFilesIn DPP-0_clean.fq --sjdbGTFfile mm39.gtf --sjdbOverhang 149 --outFileNamePrefix DPP-0-mm39 --outSAMtype BAM SortedByCoordinate
picard
Picard由下列两样构成：基于Java的来处理SAM文件的命令行工具，和一个Java API (HTSJDK) （Java应用程序接口），这个API可以用于创建可以读写SAM文件的新软件。它支持SAM文本格式和SAM二进制格式 (BAM)。
注：参考基因组序列需要.dict文件和.fai文件，可参考https://software.broadinstitute.org/gatk/documentation/article?id=1601
PS: picard已更新语法！！！（https://gatk.broadinstitute.org/hc/en-us/articles/360037052812-MarkDuplicates-Picard-）
/home/u20111230014/workspace/genome/GRCm39/STAR_GRCm39/picard_AddOrReplaceReadGroups.slurm
操作路径：/home/u20111230014/workspace/genome/GRCm39/STAR_GRCm39/
示例
java -jar ~/biosoft/picard/picard.jar CreateSequenceDictionary R=GRCm38.p5.genome.fa O=GRCm38.p5.genome.dict
samtools faidx GRCm38.p5.genome.fa
实例
操作路径：/home/u20111230014/workspace/genome/GRCm39/STAR_GRCm39
java -jar /home/public/software/miniconda2/share/picard-2.26.2-0/picard.jar CreateSequenceDictionary -R GRCm39.fa -O GRCm39.dict
samtools faidx GRCm39.fa
java -jar /home/public/software/miniconda2/share/picard-2.26.2-0/picard.jar CreateSequenceDictionary -R mm39.fa -O mm39.dict
samtools faidx mm39.fa
###
实例
java -jar /home/public/software/miniconda2/share/picard-2.26.2-0/picard.jar CreateSequenceDictionary -R mm10.fa -O mm10.dict
###比对结果文件预处理
在用STAR的2-pass mode比对时，由于考虑到后续还要给bam文件添加RG标签（GATK要求的），所以就没有在比对输出时就对其先排序，反正用picard在添加RG标签时也能顺便排序:
使用picard在添加RG标签时顺便排序
PS:GATK要求输入的bam文件包含Read groups，如果没有就会报错。
Read group是@RG开始。
示例查找Read groups
samtools view -H sample.bam | grep '@RG'
[u20111230014@cpu15 STAR_GRCm39]$ samtools view -H DPP-0-GRCm39Aligned_sortedByCoord_sorted.bam |grep '@RG'
@RG ID:DPP-0-GRCm39 LB:mRNA PL:illumina SM:DPP-0 PU:HiSeq2500
PS: 1.Picard是通过使用HTSJDK Java 库来实现的，java -jar和软件的具体路径一定要写上，
2.RGLB
JeremyL
--RGID，即-ID = Read group identifier
--RGLB，即LB= DNA preparation library identifier: 对一个read group的reads进行重复序列标记时，需要使用LB来区分reads来自那条lane;有时候，同一个库可能在不同的lane上完成测序;为了加以区分，同一个或不同库只要是在不同的lane产生的reads都要单独给一个ID。
--RGPL，即PL= Platform/technology used to produce the read, 测序使用的平台
--RGPU，即PU= Platform Unit，由三部分组成： {FLOWCELL_BARCODE}.{LANE}.{SAMPLE_BARCODE}。GATK 使用时，PU不是必须要求的；但是PU与ID同时存在时，PU优先级高于ID。
--RGSM，即SM= Sample，GATK产生的VCF文件也使用这个名字
实例
picard=/home/public/software/miniconda2/share/picard-2.26.2-0/picard.jar
java -jar $picard AddOrReplaceReadGroups -I DPP-0-GRCm39Aligned.sortedByCoord.out.bam -O DPP-0-GRCm39Aligned_sortedByCoord_sorted.bam -SO coordinate --RGID DPP-0-GRCm39 --RGLB DPP-0-GRCm39 --RGPL illumina --RGPU HiSeq2500 --RGSM DPP-0
java -jar $picard AddOrReplaceReadGroups -I DPP-1-GRCm39Aligned.sortedByCoord.out.bam -O DPP-1-GRCm39Aligned_sortedByCoord_sorted.bam -SO coordinate --RGID DPP-1-GRCm39 --RGLB DPP-1-GRCm39 --RGPL illumina --RGPU HiSeq2500 --RGSM DPP-1
##用picard套件MarkDuplicates命令对duplicate reads进行标记，并构建index （--CREATE_INDEX true ）
picard=/home/public/software/miniconda2/share/picard-2.26.2-0/picard.jar
java -jar $picard MarkDuplicates -I DPP-0-GRCm39Aligned_sortedByCoord_sorted.bam -O DPP-0-GRCm39_marked_duplicates.bam -M DPP-0-GRCm39_marked_dup_metrics.txt
--CREATE_INDEX true --VALIDATION_STRINGENCY SILENT
java -jar $picard MarkDuplicates -I DPP-1-GRCm39Aligned_sortedByCoord_sorted.bam -O DPP-1-GRCm39_marked_duplicates.bam -M DPP-1-GRCm39_marked_dup_metrics.txt
--CREATE_INDEX true --VALIDATION_STRINGENCY SILENT
参考：（https://gatk.broadinstitute.org/hc/en-us/articles/360057439771-MarkDuplicates-Picard-）
https://cncbi.github.io/Picard-Manual-CN/
[HTSJDK](https://links.jianshu.com/go?to=http%3A%2F%2Fgithub.com%2Fsamtools%2Fhtsjdk)
