31

用samtools idxstats来对de novo的转录组数据计算表达量

de novo的转录组数据,比对的时候一般用的是自己组装好的trinity.fasta序列(挑选最长蛋白的转录本序列)来做参考,用bowtie2等工具直接将原始序列比对即可。所以比对 sam/bam文件本身就包含了参考序列的每一条转录本序列ID,直接对 sam/bam文件进行counts就知道每一个基因的表达量啦!

本来我是准备自己写脚本对sam文件进行counts就好,但是发现了samtools自带这样的工具:http://www.htslib.org/doc/samtools.html 

如果是针对基因组序列,那么这个功能用处不大,但是针对转录本序列,统计出来的就是我们想要的转录本表达量。 Continue reading

21

自学无参RNAseq数据分析第一讲之参考文献解读

这是我为新创办的 生信技能树 论坛写的帖子,也适合本博客,所以转载过来: http://www.biotrainee.com/thread-243-1-1.html 

以前做的都是有参转录组分析,只需要找到参考基因组和注释文件,然后走QC-->alignment-->counts->DEG-->annotation的流程图即可。
现在开始学习新的东西了,就是无参转录组分析,这里记录一下自己的学习笔记,首先还是资料收集,这次,我就针对性的看5个 全流程化的转录组 de novo 分析 文章,如下:
http://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-15-554  2014年栀子花的花瓣衰老的标准de novo 转录组分析,数据如下:用Trinity做组装,用NCBI non-redundant (Nr) database库做注释,做了差异分析(栀子花花期分成4个阶段),GO/KEGG注释,然后做了RT-qPCR的实验验证。
多做了一个 Clusters of Orthologus Groups (COG)的数据库注释

Raw Reads
Clean Reads
Contigs
Unigenes
Annotated
Transcriptome
55,092,396
50,335,672
102,263
57,503
39,459

 

http://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-15-236  2014 巴西橡胶树的研究,是一个综合多组织样本的RNA库,ployT建库,454测序,用的是est2Assembly 和gsassembler 软件做组装,用 NCBI RefSeq, Plant Protein Database 做注释,因为没有分组,所以不必做差异分析,只需要找SNV和SSR标记即可,最后也是做GO/KEGG注释

https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-016-2633-2 2015 萝卜,用illumina进行转录组测序,用Trinity组装,用RPKM值算unigene的表达量,也是用 BLASTx来对Trinity结果进行注释,注释到NR,NT,Swiss-Prot,GO,COG,kegg数据库,其中GO注释用的是Blast2GO,最后也做了RT-qPCR 实验验证,某些基因在leaf里面的表达量显著高于其它tissue,有原始数据:http://www.ncbi.nlm.nih.gov/sra/?term=SRX1671013
转录组分析结果结果:A total of 54.64 million clean reads and 111,167 contigs representing 53,642 unigenes were obtained from the radish leaf transcriptome.

http://www.nature.com/articles/srep08259 2015 芹菜 叶片发育中木质素的探究,测序的reads是A total of 32,477,416 quality reads were recorded for the leaves at Stage 1, 53,675,555 at Stage 2, and 27,158,566 at Stage 3, respectively.,也是用Trinity组装,kmer值设为25,组装结果:33,213 unigenes with an average length of 1,478 bp, a maximum length of 17,075 bp, and an N50 of 2,060 bp,然后用eggNOG/GO/KEGG数据库来注释。文章正文给了所用到的软件和数据库的详细链接
最后还用了 real-time PCR assays          来看 roots, stems, petioles, and leaf blade 这些组织的基因表达差异情况

http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0128659 对 三疣梭子蟹 的卵巢和睾丸的转录组研究,,也是标准的转录组de novo 分析流程,非常值得借鉴
NCBI有上传原始数据:SRR1920180  和SRR1920180  

总结好这5篇文献的数据分析流程,就差不多明白如何做无参的转录组de novo分析了

十一 01

WES(七)看de novo变异情况

de novo变异寻找这个也属于snp-calling的一部分,但是有点不同的就是该软件考虑了一家三口的测序文件,找de novo突变

软件有介绍这个功能:http://varscan.sourceforge.net/trio-calling-de-novo-mutations.html

而且还专门有一篇文章讲ASD和autism与de novo变异的关系,但是文章不清不楚的,没什么意思

Trio Calling for de novo Mutations

image001

Min coverage:   10

Min reads2:     4

Min var freq:   0.2

Min avg qual:   15

P-value thresh: 0.05

Adj. min reads2:        2

Adj. var freq:  0.05

Adj. p-value:   0.15

Reading input from trio.filter.mpileup

1371416525 bases in pileup file (137M的序列)

83123183 met the coverage requirement of 10 (其中有83M的测序深度大于10X)

145104 variant positions (132268 SNP, 12836 indel) (共发现15.5万的变异位点)

4403 were failed by the strand-filter

139153 variant positions reported (126762 SNP, 12391 indel)

502 de novo mutations reported (376 SNP, 126 indel) (真正属于 de novo mutations只有502个)

1734 initially DeNovo were re-called Germline

12 initially DeNovo were re-called MIE

3 initially DeNovo were re-called MultAlleles

522 initially MIE were re-called Germline

1 initially MIE were re-called MultAlleles

3851 initially Untransmitted were re-called Germline

然后我看了看输出的文件trio.mpileup.output.snp.vcf

软件是这样解释的:The output of the trio subcommand is a single VCF in which all variants are classified as germline (transmitted or untransmitted), de novo, or MIE.

  • FILTER - mendelError if MIE, otherwise PASS
  • STATUS - 1=untransmitted, 2=transmitted, 3=denovo, 4=MIE
  • DENOVO - if present, indicates a high-confidence de novo mutation call

里面的信息量好还是不清楚

我首先对我们拿到的trio.de_novo.mutaion.snp.vcf文件进行简化,只看基因型!
head status.txt   (顺序是dad,mom,child)
STATUS=2 0/0 0/1 0/1
STATUS=2 1/1 1/1 1/1
STATUS=2 0/1 0/0 0/1
STATUS=2 1/1 1/1 1/1
STATUS=1 0/1 0/0 0/0
STATUS=1 0/1 0/0 0/0
STATUS=2 1/1 1/1 1/1
STATUS=2 1/1 1/1 1/1
STATUS=2 1/1 1/1 1/1
STATUS=2 0/1 0/1 0/1
那么总结如下:
  26564 STATUS=1 无所以无 (0/0 0/1 0/0或者 0/1 0/0 0/0等等)
  97764 STATUS=2 有所以有 (1/1 1/1 1/1 或者0/1 0/1 1/1等等)
    385 STATUS=3 无中生有 (0/0 0/0 0/1 或者0/0 0/0 1/1)
   1485 STATUS=4 有中生无 (1/1 0/1 0/0 等等)

我用annovar注释了一下

/home/jmzeng/bio-soft/annovar/convert2annovar.pl -format vcf4  trio.mpileup.output.snp.vcf > trio.snp.annovar

/home/jmzeng/bio-soft/annovar/annotate_variation.pl -buildver hg19  --geneanno --outfile  trio.snp.anno trio.snp.annovar /home/jmzeng/bio-soft/annovar/humandb

结果是:

A total of 132268 locus in VCF file passed QC threshold, representing 132809 SNPs (86633 transitions and 46176 transversions) and 3 indels/substitutions

可以看到最后被注释到外显子上面的突变有两万多个

23794  284345 3123333 trio.snp.anno.exonic_variant_function

这个应该是非常有意义的,但是我还没学会后面的分析。只能先做到这里了