GATK入门

1. Add read groups (Picard tools)
   AddOrReplaceReadGroups.jar I=sorted.bam_file O=s1.rg.bam RGLB=genome RGPL=ILLUMINA
   RGPU=GATKv4 RGSM=sample_name VALIDATION_STRINGENCY=LENIENT
   用了picard tools的AddOrReplaceReadGroups.jar加了read group，根据以前百度GATK的经验，了解GATK要求bam文件的header必须包含@RG，所以这一步应该是前面比对时候，没有在参数中增加相应部分。所以如果我在比对的时候增加了这个参数，这一步就可以免了。
   
   
   
2. Mark duplicates (Picard tools)
   MarkDuplicates.jar INPUT=s1.rg.bam OUTPUT=s2.dedup.bam ASSUME_SORTED=TRUE
   VALIDATION_STRINGENCY=LENIENT METRICS_FILE=s2.dedup.metrics
   
   

用picard tools的MarkDuplicates.jar去重。这是因为二代测序有一部桥式PCR扩增底物的过程，在100x以上出现reads重复，很大可能是PCR扩增重复了，所以可以直接去掉。

3. Index (samtools)
   samtools index s2.dedup.bam
   samtools的index，建立索引，提高检索效率。
   
   
4. Realign reads (create intervals first, then do IndelRealigner) (GATK)
   GenomeAnalysisTK.jar -I s2.dedup.bam -R ref_file -T RealignerTargetCreator -o s3.intervals
   GenomeAnalysisTK.jar -T IndelRealigner -I s2.dedup.bam -R ref_file -targetIntervals s3.intervals -o
   s4.realn.bam
   先建立indel区间文件，然后对该区域进行reads重比对。这一步我没有经验，然后我在GATK的论坛搜索了这个工具，发现原因是indel会导致附近的错配，所以需要借由这一步降低indel附近的假阳性。但是最新的HaplotypeCaller or MuTect2已经不需要了，这些工具的haplotype组装步奏效果类似。这里的UnifiedGenotyper or the original MuTect.还是要的。
   
   
   
   
5. Unified genotyper (GATK)
   GenomeAnalysisTK.jar -T UnifiedGenotyper -R ref_file -I s4.realn.bam -glm BOTH -o s5.UG1.vcf -mbq 30 -nt
   使用UnifiedGenotyper寻找variant。论坛搜索发现，现在推荐用HaplotypeCaller，效果更好。
   4
   
   
   
6. Base score recalibrator (GATK)
   GenomeAnalysisTK.jar -T BaseRecalibrator -I s4.realn.bam -R ref_file -knownSites s5.UG1.vcf -o s6.recal
   

根据上一步得到vcf文件对第四步得到BAM文件进行碱基重校准，得到校准所需的文件。

7. Print Reads (GATK)
   GenomeAnalysisTK.jar -T PrintReads -R ref_file -I s4.realn.bam -BQSR s6.recal -o s7.recal.bam
   
8. Unified Genotype (GATK)
   GenomeAnalysisTK.jar -T UnifiedGenotyper -R ref_file -I s7.recal.bam -glm BOTH -o s8.UG2.vcf -mbq 30
   
9. Base score recalibrator (GATK)
   GenomeAnalysisTK.jar -T BaseRecalibrator -I s4.realn.bam -R ref_file -knownSites s8.UG1.vcf -o s9.recal
   
10. Print Reads (GATK)
    GenomeAnalysisTK.jar -T PrintReads -R ref_file -I s4.realn.bam -BQSR s9.recal -o s10.recal.bam
    
11. Unified Genotyper (GATK)
    GenomeAnalysisTK.jar -T UnifiedGenotyper -R ref_file -I s10.recal.bam -glm BOTH
    

基本方法就是以生物学意义的方式计算基因表达量，然后通过统计学分析表达量寻找具有统计学显著性差异的基因，从而

说人话就是，基因A和基因B的平均值之差与两者中较小的比值。选择2-3倍的基因作为结果（为什么是2-3倍，就是大家约定俗成）。

用到的数据来自于一篇文献“A pooled shRNA screen for regulators of primary mammary stem and progenitor cells identifies roles for Asap1 and Prox1”。

install.packages(“R.utils”)

read.delim(files[1], nrow=5)

EntrezID GeneLength Count

sample1 <- read.delim(files[1],header = T)[,c(1,3)]

count.table <- data.frame(sample1)

然后要提供colData,其中colData存放列名要和countData的行名相同。

colData

all(rownames(colData) %in% colnames(count.table))

最后导入数据:

到底用啥：数据集小于30 -> rlog，大数据集 -> VST。还有这个处理过程不是用于差异检验的，在DESeq分析中会自动选择最合适的所以你更不需要纠结了，记得用raw count。

dds <- estimateSizeFactors(dds)

colnames(df)[1:2] <- c(“x”, “y”)

结果就是转换后更加集中了。

同样的可以用Poisson Distance (Witten 2011)计算距离，计算方式如下：

plotPCA(rld,intgroup=c(“condition”))

能够明显的发现不同处理的距离离得很远。

比如可以指定比较对象Basal和PL,可以用mcols查看结果存储的元数据，了解列名的含义。

于是乎结果就和limma差不多了。

为什么要做多重试验矫正？

如果我们选择了所有小于或等于矫正p-value阈值的所有显著性基因，假阳性比例（ false discovery rate, FDR）是多少？

Up

head(resSig[ order(resSig$log2FoldChange), ])

Down

head(resSig[ order(resSig$log2FoldChange, decreasing = TRUE), ])

heatmap

大致可以发现同一组的基因颜色是相同的，也就是说表达量相近。

首先对基因表达分析（上）做一个简单的回顾

然后用Salmon建立索引：

salmon index -t athal.fa.gz -i athal_index

tximport可以纠正不同样本基因长度的潜在改变（比如说differential isoform usage）；能够用于导入 (Salmon, Sailfish, kallisto)程序的输出文件；能够避免丢弃那些比对到多个基因的同源序列，从而提高灵敏度

install.packages(“readr”)

install.packages(“rsjon”)

reading in files with read_tsv

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

summarizing abundance

summarizing counts

summarizing length

head(txi$lenght)head(txi$counts)

我希望通过循环的方式，以每个matrix的两列作为输入做线性拟合, 如下

lm(u[,1]~u[,2])

我忘记了c()会把所有matrix降成一维，这就很尴尬了。这是因为R不支持对非向量集合的循环操作。如果想实现非向量集合的循环，只能曲线救国，如lapply或get，这里介绍get。

for (m in c(“u”,“v”)){

• z <- get(m) 

• print(z) 

• }
  [,1] [,2]
  [1,] 1 9
  [2,] 2 10
  [3,] 3 11
  [4,] 4 12
  [5,] 5 13
  [6,] 6 14
  [7,] 7 15
  [8,] 8 16
  [,1] [,2]
  [1,] 17 23
  [2,] 18 24
  [3,] 19 25
  [4,] 20 26
  [5,] 21 27
  [6,] 22 28
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  

于是，这解决了我昨天晚上遇到的问题。我打算对3个DESeq2的比较table进行过滤，

本来代码长下面这个样子，各种他提示错误：

所以最后的代码长这个样子：

其实可以用list，然后用lapply, 然后分别提取的。

resSigUp <- lapply(list(res1,res2,res3), subset, padj < 0.01 & log2FoldChange >0 )