转录组分析---Hisat2+StringTie+Ballgown使用

时间:2021-02-17 21:08:32

转录组分析---Hisat2+StringTie+Ballgown使用

(2016-10-10 08:14:45)

标签:

生物信息学

转录组

 
1.Hisat2建立基因组索引:
First, using the python scripts included in the HISAT2 package, extract splice-site and exon information from the gene
annotation file:
 
$ extract_splice_sites.py gemome.gtf >genome.ss#得到剪接位点信息
$ extract_exons.py genome.gtf >genome.exon#得到外显子信息
 
Second, build a HISAT2 index:
 
$ hisat2-build --ss genome.ss --exon genome.exon genome.fa genome
 
备注extract_splice_sites.py 和 extract_exons.py 在hisat2软件包中涵盖了,这两步不是必须的,只是为了发现剪切位点,也可以直接:
$ hisat2-build  genome.fa genome
 
2. 利用hisat2比对到基因组:
 
hisat2 -p 8 --dta -x genome -1 file1_1.fastq.gz -2 file1_2.fastq.gz -S file1.sam
hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1 file2_1.fastq.gz -2 file2_2.fastq.gz -S file2.sam
 
备注:--dta:输出转录组型的报告文件
-x:基因组索引
-S : 输出sam文件
-p: 线程数
其他参数:
Input:
  -q                 query input files are FASTQ .fq/.fastq (default)
  --qseq             query input files are in Illumina's qseq format
  -f                 query input files are (multi-)FASTA .fa/.mfa
  -r                 query input files are raw one-sequence-per-line
  -c                 , , are sequences themselves, not files
  -s/--skip    skip the first reads/pairs in the input (none)
  -u/--upto    stop after first reads/pairs (no limit)
  -5/--trim5   trim bases from 5'/left end of reads (0)
  -3/--trim3   trim bases from 3'/right end of reads (0)
  --phred33          qualities are Phred+33 (default)
  --phred64          qualities are Phred+64
  --int-quals        qualities encoded as space-delimited integers
 
 Alignment:
  -N           max # mismatches in seed alignment; can be 0 or 1 (0)
  -L           length of seed substrings; must be >3, <32 (22)
  -i          interval between seed substrings w/r/t read len (S,1,1.15)
  --n-ceil    func for max # non-A/C/G/Ts permitted in aln (L,0,0.15)
  --dpad       include extra ref chars on sides of DP table (15)
  --gbar       disallow gaps within nucs of read extremes (4)
  --ignore-quals     treat all quality values as 30 on Phred scale (off)
  --nofw             do not align forward (original) version of read (off)
  --norc             do not align reverse-complement version of read (off)
 
 Spliced Alignment:
  --pen-cansplice              penalty for a canonical splice site (0)
  --pen-noncansplice           penalty for a non-canonical splice site (12)
  --pen-canintronlen          penalty for long introns (G,-8,1) with canonical splice sites
  --pen-noncanintronlen       penalty for long introns (G,-8,1) with noncanonical splice sites
  --min-intronlen              minimum intron length (20)
  --max-intronlen              maximum intron length (500000)
  --known-splicesite-infile   provide a list of known splice sites
  --novel-splicesite-outfile  report a list of splice sites
  --novel-splicesite-infile   provide a list of novel splice sites
  --no-temp-splicesite               disable the use of splice sites found
  --no-spliced-alignment             disable spliced alignment
  --rna-strandness          Specify strand-specific information (unstranded)
  --tmo                              Reports only those alignments within known transcriptome
  --dta                              Reports alignments tailored for transcript assemblers
  --dta-cufflinks                    Reports alignments tailored specifically for cufflinks
 
 Scoring:
  --ma         match bonus (0 for --end-to-end, 2 for --local)
  --mp ,   max and min penalties for mismatch; lower qual = lower penalty <2,6>
  --sp ,   max and min penalties for soft-clipping; lower qual = lower penalty <1,2>
  --np         penalty for non-A/C/G/Ts in read/ref (1)
  --rdg ,  read gap open, extend penalties (5,3)
  --rfg ,  reference gap open, extend penalties (5,3)
  --score-min min acceptable alignment score w/r/t read length
                     (L,0.0,-0.2)
 
 Reporting:
  (default)          look for multiple alignments, report best, with MAPQ
   OR
  -k           report up to alns per read; MAPQ not meaningful
   OR
  -a/--all           report all alignments; very slow, MAPQ not meaningful
 
 Effort:
  -D           give up extending after failed extends in a row (15)
  -R           for reads w/ repetitive seeds, try sets of seeds (2)
 
 Paired-end:
  --fr/--rf/--ff     -1, -2 mates align fw/rev, rev/fw, fw/fw (--fr)
  --no-mixed         suppress unpaired alignments for paired reads
  --no-discordant    suppress discordant alignments for paired reads
 
 Output:
  -t/--time          print wall-clock time taken by search phases
  --un           write unpaired reads that didn't align to
  --al           write unpaired reads that aligned at least once to
  --un-conc      write pairs that didn't align concordantly to
  --al-conc      write pairs that aligned concordantly at least once to
  (Note: for --un, --al, --un-conc, or --al-conc, add '-gz' to the option name, e.g.
  --un-gz , to gzip compress output, or add '-bz2' to bzip2 compress output.)
  --quiet            print nothing to stderr except serious errors
  --met-file  send metrics to file at (off)
  --met-stderr       send metrics to stderr (off)
  --met        report internal counters & metrics every secs (1)
  --no-head          supppress header lines, i.e. lines starting with @
  --no-sq            supppress @SQ header lines
  --rg-id     set read group id, reflected in @RG line and RG:Z: opt field
  --rg        add ("lab:value") to @RG line of SAM header.
                     Note: @RG line only printed when --rg-id is set.
  --omit-sec-seq     put '*' in SEQ and QUAL fields for secondary alignments.
 
 Performance:
  -o/--offrate override offrate of index; must be >= index's offrate
  -p/--threads number of alignment threads to launch (1)
  --reorder          force SAM output order to match order of input reads
  --mm               use memory-mapped I/O for index; many 'bowtie's can share
 
 Other:
  --qc-filter        filter out reads that are bad according to QSEQ filter
  --seed       seed for random number generator (0)
  --non-deterministic seed rand. gen. arbitrarily instead of using read attributes
  --version          print version information and quit
  -h/--help          print this usage message
 
 
3.  将sam文件sort并转化成bam:
 
$ samtools sort -@ 8 -o file1.bam file1.sam
$ samtools sort -@ 8 -o file2.bam file2.sam
 
4. 组装转录本:
 
$ stringtie -p 8 -G genome.gtf -o file1.gtf –l file1 file1.bam
$ stringtie -p 8 -G genome.gtf -o file2.gtf –l file2 file2.bam
lncRNA (-f 0.01 -a 10 -j 1 -c 0.01)
其中:
 -G reference annotation to use for guiding the assembly process (GTF/GFF3)
 -l name prefix for output transcripts (default: STRG)
 -f minimum isoform fraction (default: 0.1)
 -m minimum assembled transcript length (default: 200)
 -o output path/file name for the assembled transcripts GTF (default: stdout)
 -a minimum anchor length for junctions (default: 10)
 -j minimum junction coverage (default: 1)
 -t disable trimming of predicted transcripts based on coverage
    (default: coverage trimming is enabled)
 -c minimum reads per bp coverage to consider for transcript assembly
    (default: 2.5)
 -v verbose (log bundle processing details)
 -g gap between read mappings triggering a new bundle (default: 50)
 -C output a file with reference transcripts that are covered by reads
 -M fraction of bundle allowed to be covered by multi-hit reads (default:0.95)
 -p number of threads (CPUs) to use (default: 1)
 -A gene abundance estimation output file
 -B enable output of Ballgown table files which will be created in the
    same directory as the output GTF (requires -G, -o recommended)
 -b enable output of Ballgown table files but these files will be
    created under the directory path given as
 -e only estimate the abundance of given reference transcripts (requires -G)
 -x do not assemble any transcripts on the given reference sequence(s)
 -h print this usage message and exit
 
 
5. 合并所有样本的gtf文件
 
$ stringtie --merge -p 8 -G genome.gtf -o stringtie_merged.gtf mergelist.txt
 
6. 新转录本的注释(lncRNA必备,普通转录组忽略)
 
gffcompare –r genomegtf –G –o merged stringtie_merged.gtf
 
备注:gffcompare 是独立软件,下载地址http://ccb.jhu.edu/software/stringtie/gff.shtml,结果如下;
= Predicted transcript has exactly the same introns as the reference transcript
c Predicted transcript is contained within the reference transcript
j Predicted transcript is a potential novel isoform that shares at least one splice junction with a reference transcript
e Predicted single-exon transcript overlaps a reference exon plus at least 10 bp of a reference intron, indicating a possible pre-mRNA fragment
i Predicted transcript falls entirely within a reference intron
o Exon of predicted transcript overlaps a reference transcript
p Predicted transcript lies within 2 kb of a reference transcript (possible polymerase run-on fragment)
r Predicted transcript has >50% of its bases overlapping a soft-masked (repetitive) reference sequence
u Predicted transcript is intergenic in comparison with known reference transcripts
x Exon of predicted transcript overlaps reference but lies on the opposite strand
s Intron of predicted transcript overlaps a reference intron on the opposite strand
 
7. 转录本定量和下游ballgown软件原始文件构建:
 
$ stringtie –e –B -p 8 -G stringtie_merged.gtf -o ballgown/file1/file1.gtf file1.bam
$ stringtie –e –B -p 8 -G stringtie_merged.gtf -o ballgown/file2/file2.gtf file2.bam
 
8. Ballgown差异表达分析:
 
>library(ballgown)
>library(RSkittleBrewer)
>library(genefilter)
>library(dplyr)
>library(devtools)
>pheno_data = read.csv("geuvadis_phenodata.csv")#读取表型数据
>bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "file", pData=pheno_data)#读取表达量
>bg_chrX_filt = subset(bg_chrX,"rowVars(texpr(bg_chrX)) >1",genomesubset=TRUE)#过滤低表达量基因
>results_transcripts = stattest(bg_chrX_filt,feature="transcript",covariate="sex",adjustvars =c("population"), getFC=TRUE, meas="FPKM")#差异表达分析,运用的是一般线性模型,比较组sex,影响因素:population
>results_genes = stattest(bg_chrX_filt, feature="gene",covariate="sex", adjustvars = c("population"), getFC=TRUE,meas="FPKM")#基因差异表达
>results_transcripts=data.frame(geneNames=ballgown::geneNames(bg_chrX_filt),geneIDs=ballgown::geneIDs(bg_chrX_filt), results_transcripts)#增加基因名字,id
>results_transcripts = arrange(results_transcripts,pval)#按pval sort
>results_genes = arrange(results_genes,pval)
>write.csv(results_transcripts, "chrX_transcript_results.csv",
row.names=FALSE)
>write.csv(results_genes, "chrX_gene_results.csv",
row.names=FALSE)
>subset(results_transcripts,results_transcripts$qval<0.05)
>subset(results_genes,results_genes$qval<0.05)
 
9. 结果可视化:
 
>tropical= c('darkorange', 'dodgerblue',
'hotpink', 'limegreen', 'yellow')
>palette(tropical)
>fpkm = texpr(bg_chrX,meas="FPKM")
>fpkm = log2(fpkm+1)
>boxplot(fpkm,col=as.numeric(pheno_data$sex),las=2,ylab='log2(FPKM+1)')
>ballgown::transcriptNames(bg_chrX)[12]
## 12
## "NM_012227"
>ballgown::geneNames(bg_chrX)[12]
## 12
## "GTPBP6"
>plot(fpkm[12,] ~ pheno_data$sex, border=c(1,2),
main=paste(ballgown::geneNames(bg_chrX)[12],' : ',
ballgown::transcriptNames(bg_chrX)[12]),pch=19, xlab="Sex",
ylab='log2(FPKM+1)')
>points(fpkm[12,] ~ jitter(as.numeric(pheno_data$sex)),
col=as.numeric(pheno_data$sex))
>plotTranscripts(ballgown::geneIDs(bg_chrX)[1729], bg_chrX, main=c('Gene XIST in sample ERR188234'), sample=c('ERR188234'))
>plotMeans('MSTRG.56', bg_chrX_filt,groupvar="sex",legend=FALSE)