!!!對(duì) RNA-seq 產(chǎn)出的數(shù)據(jù)進(jìn)行變異檢測(cè)分析，與常規(guī)重測(cè)序的主要區(qū)別就在序列比對(duì)這一步，因?yàn)?RNA-seq 的數(shù)據(jù)是來自轉(zhuǎn)錄本的，比對(duì)到參考基因組需要跨越轉(zhuǎn)錄剪切位點(diǎn)，所以 RNA-seq 進(jìn)行變異檢測(cè)的重點(diǎn)就在于跨剪切位點(diǎn)的精確序列比對(duì)。
STAR相比較其他兩款軟件有較高的唯一比對(duì)率；STAR會(huì)將沒有paired mapping上的reads都剔除，避免single reads比對(duì)到基因組上；并且STAR對(duì)lower-quality（包括more soft-clipped和錯(cuò)配堿基）比對(duì)有較高的容忍度
構(gòu)建基因組索引
STAR --runThreadN 25 --runMode genomeGenerate --genomeDir index_dir --genomeFastaFiles genome.fasta --sjdbGTFfile genome.gtf --sjdbOverhang 149
--runThreadN：線程數(shù)。
--runMode genomeGenerate：構(gòu)建基因組索引。
--genomeDir：索引目錄。（index_dir一定要是存在的文件夾，需提前建好）
--genomeFastaFiles：基因組文件。
--sjdbGTFfile：基因組注釋文件。
--sjdbOverhang：reads長(zhǎng)度減1。
實(shí)例
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
#若出現(xiàn)報(bào)錯(cuò)如下：
[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
##解決
一般平臺(tái)使用slurm時(shí)，cpu核數(shù)與內(nèi)存之間要相對(duì)平衡，建議1:5，即#SBATCH -c 2; #SBATCH --mem=10G。
將腳本的cpu改大一點(diǎn)，并配上 --limitGenomeGenerateRAM 報(bào)錯(cuò)的內(nèi)存數(shù)字
#!/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
為了發(fā)現(xiàn)更加靈敏的new junction，STAR建議使用2-pass mode，其能增加檢測(cè)到的new junction數(shù)目，使得更多的splices reads能mapping到new junction。因此STAR先用一般參數(shù)做一遍mapping，收集檢測(cè)到的junction信息，然后利用這已經(jīng)annotated junction來做第二次mapping
雙端比對(duì)
示例
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
##單端比對(duì)
示例
STAR \
--runThreadN 20 \
--genomeDir hg19_STAR_db \
--readFilesIn reads.fq \
--sjdbGTFfile hg19.gtf \
--sjdbOverhang 100 \
--outFileNamePrefix sampleA \
--outSAMtype BAM SortedByCoordinate
##命令參數(shù)也很好理解：
--runThreadN ：設(shè)置線程數(shù)
--genomeDir ：index輸出的路徑
--genomeFastaFiles ：參考基因組序列
--sjdbGTFfile ：參考基因組注釋文件
--sjdbOverhang ：這個(gè)是reads長(zhǎng)度的最大值減1，默認(rèn)是100
如果是per-sample(單樣本)的2-pass mapping，可以直接用--twopassMode Basic參數(shù)將第兩步mapping中的make index省去，直接再mapping
操作路徑：/home/u20111230014/workspace/genome/STAR_mm10
實(shí)例
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
實(shí)例
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
實(shí)例
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由下列兩樣構(gòu)成：基于Java的來處理SAM文件的命令行工具，和一個(gè)Java API (HTSJDK) （Java應(yīng)用程序接口），這個(gè)API可以用于創(chuàng)建可以讀寫SAM文件的新軟件。它支持SAM文本格式和SAM二進(jìn)制格式 (BAM)。
注：參考基因組序列需要.dict文件和.fai文件，可參考https://software.broadinstitute.org/gatk/documentation/article?id=1601
PS: picard已更新語法?。。。?a target="_blank">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
實(shí)例
操作路徑：/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
###
實(shí)例
java -jar /home/public/software/miniconda2/share/picard-2.26.2-0/picard.jar CreateSequenceDictionary -R mm10.fa -O mm10.dict
###比對(duì)結(jié)果文件預(yù)處理
在用STAR的2-pass mode比對(duì)時(shí)，由于考慮到后續(xù)還要給bam文件添加RG標(biāo)簽（GATK要求的），所以就沒有在比對(duì)輸出時(shí)就對(duì)其先排序，反正用picard在添加RG標(biāo)簽時(shí)也能順便排序:
使用picard在添加RG標(biāo)簽時(shí)順便排序
PS:GATK要求輸入的bam文件包含Read groups，如果沒有就會(huì)報(bào)錯(cuò)。
?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 庫來實(shí)現(xiàn)的，java -jar和軟件的具體路徑一定要寫上，
2.RGLB
JeremyL
--RGID，即-ID = Read group identifier
--RGLB，即LB= DNA preparation library identifier: 對(duì)一個(gè)read group的reads進(jìn)行重復(fù)序列標(biāo)記時(shí)，需要使用LB來區(qū)分reads來自那條lane;有時(shí)候，同一個(gè)庫可能在不同的lane上完成測(cè)序;為了加以區(qū)分，同一個(gè)或不同庫只要是在不同的lane產(chǎn)生的reads都要單獨(dú)給一個(gè)ID。
--RGPL，即PL= Platform/technology used to produce the read, 測(cè)序使用的平臺(tái)
--RGPU，即PU= Platform Unit，由三部分組成： {FLOWCELL_BARCODE}.{LANE}.{SAMPLE_BARCODE}。GATK 使用時(shí)，PU不是必須要求的；但是PU與ID同時(shí)存在時(shí)，PU優(yōu)先級(jí)高于ID。
--RGSM，即SM= Sample，GATK產(chǎn)生的VCF文件也使用這個(gè)名字
實(shí)例
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命令對(duì)duplicate reads進(jìn)行標(biāo)記，并構(gòu)建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)
