2017.3.24 上午

==================== 基因组 2017.3.24

1.什么样的物种适合全基因组测序？能被高引用的，才可能高分文章。天山雪莲，从稀有、中药角度，老外估计不喜欢。如果有天山雪莲的体外培养系统，做成耐寒的模式生物，则有可能发高分。 end9:45

2.技术流、思路哪个重要？几k的测序结果，就一个进化树，为什么发到CNS？还没有自己的真菌分析复杂、透彻？因为自己的真菌没人关注，而别人的传染病病毒则很多人关注。

比如SARS时，很多机构都很关注：知名实验室、制药公司、政府部门。华大当时通过私人关系“偷”的菌株，测序的。

影响力由什么决定？模式生物、病原性、经济性状、进化地位。决定于有多少读者群。

end9:51

3.课题设计很重要？只有基因组可能不够，是静态的。多层次：转录、表观、代谢。只要是动态的，就涉及取样问题，要提前设计。

经济方向怎么解释？以前一直做雪莲药用的，次生代谢基础做得好。但是读者可能不关注：做基因组测序前，同行调研。抗逆，可能关注度排第一。 end10:00

4.基因组测序：二代+三代当前是主流。很多物种测序过了，怎么发掘新的模式生物？已测序物种：https://www.ncbi.nlm.nih.gov/guide/genomes-maps/ 当该生物在进化树分支上、小分支上首发时，paper至少一半要分析进化。进化的魅力：所有结论都是推论。无法证伪。

如果有材料、没idea，销售一般推销什么？自己公司擅长的，而不是该材料最适合的。所以设计很重要。 end10:14

5.组学文章3大块：测序、关键遗传性状、比较基因组学分析。

1)half de novo 测序：近亲物种测序过了，该物种可能只有SNP、SV等少量差异。 2)重测序适用于：该物种测序过了， 2013橡胶基因组很碎，剩8000个scaffold。再做时，不能不管第一个报道，剩下2000个scaffold。 3)population测序：第一株高精度de novo，其他做resequnce，和第一株比对。

组学文章最忌讳：一篇文章什么都讲 = 什么都没讲。《我的一天》作文不能使流水账，必须详细写一件事。要讲最吸引人的部分，1-2个经济性状。结合转录组、代谢组，动态分析。韩国辣椒基因组，转录组中辣与不辣落到一个基因上。 end10:20

6.全基因组SSR都找出来了。意味着，后面育种都要引用了。 SSR(Simple Sequence Repeats)标记是近年来发展起来的一种以特异引物PCR为基础的分子标记技术，也称为微卫星DNA(MicrosatelliteDNA),是一类由几个核苷酸（一般为1~6个）为重复单位组成的长达几十个核苷酸的串联重复序列。由于每个SSR两侧的序列一般是相对保守的单拷贝序列。

1）如果基因组碎片多，别做大进化分析，别和其他物种做比对。共线性不等于同源。碎片多，从生物学意义入手：xx基因重要，xx基因家族重要。

2）如果拼接的好，要和水稻12条染色体比对，拟南芥比对，凸显全基因组拼接的好。基因家族拷贝、收缩。

测序低价都是有原因的。 end10:35

7.怎么学生信？跟着一个完整项目。要及时问公司要文档、参数，防止人员离职。

end10:45

==================== 转录组 2017.3.24

8.转录组是必须做！有取样问题，安排时间。 de novo要三套转录组：

多组织混合样品 => 获取全长cDNA。最可靠的是蛋白组学数据，但蛋白组水平低。用于：辅助拼接、做预测用训练集。技术：长读长为主pacBio，几k，几十k，不需要打断，准确性比拼接好。三代缺点: 通量小，低丰度检测不到。二代也要获取长度长：300PE测序MySeq，先拼接成500bp片段，然后map。
多组织（根茎叶花果实、心肝脾肺肾）【深度】RNAseq表达谱 => 组织特异性表达。真核1-2G太小，推荐4-6G；
要有目的的多点取样。生物学亮点：两点取样回答问题有限，只能分析上下调基因。

end11:13

9.转录组需要精心设计。调控网络：表观。

end11:16

10.平台选择：3+2长短结合方式，仅二代也要长短结合。细菌：必须上三代，而用二代很难成环。补漏洞时间很长，一年-四年！三代没有PCR，没有GC偏差。缺点是，错误率高polyA、polyT。用二代补救。真菌：同上。

三代pacBio为拼接contig，40-60x测序。
最忌讳是一次全出来数据，要出一部分数据立刻分析，不符合预期立刻改方案。

反例：1000万，7年，基因组还没拼接好，很多漏洞，可能就是公司没水平。

end11:27

11.怎么制定基因组测序方案？ 1）花2-3万先测基因组大小？？？Genome Survey. 20x-50xPE，200bp-400bp的小库，不要500bp的大库。 jellyfish软件可以画一个kmer分布图。 kmer是统计概念。意味着：100bp read有多少17kmer？1-17,2-18,3-19(100-17+1个kmer) x：某个kmer出现的次数。 y：评估genome survey的方法。峰值/ =基因组大小。

1.2）昆虫、藻类一般都是杂合子。种群越大、寿命越短，则杂合度越多。

杂合度高的短reads，很难拼接！比重复序列还难。杂合的，数据量上升，400bp contig不会再长了。多倍体也很难。例：六倍体小麦异源杂合。

放大镜才能看到的昆虫，口腔、肠道、体表不可避免基因组污染。测序结果和所有已知基因组（植物、微生物）比对，去掉后再看kmer图。如果还有低频区域杂峰，说明还有基因组污染（未知微生物基因组）。

主峰：不是单峰，是杂合。 x轴很长，说明基因组内有很多重复序列。 x轴很长，y坐标又很大，则说明更难。

苹果的杂合度太高，只能选相对杂合度低的。虽然不是种植最广泛的品种。 kmer没有峰怎么办？测序量不够。

2）物理图谱：BioNano 等确定序列。做好数据指控。QC对后续假阳、假阴影响很大。一般建议，精准医学要用最严格的QC，宁可丢一半序列，也要只接收最高质量的数据。基因组拼接：(Genome res.2010 Feb;20(s):265-72)

形成contig：连续、无空洞
拼接scaffold：中间是有洞的
补洞：gap filling/ gap closure. 软件：
ABySS
ALLPATHS-LG
Platanus(拼接高杂合)：内部高压缩成一条序列，
SOAPdenovo

3）如果有种植资源，可以做genotyping by sequencing,GBS,构建遗传连锁群。不建议选遗传关系特别近的亲本。 4）确定连锁群在一号染色体，只能用FISH定位到染色体。

2017.3.24 pm

12. 如何评价拼接质量？

N50/N90，拼接到一半时scafold的长度。把contig 和 scaffold 从长到短进行排列,然后相加,当恰好加到1M的50％,也就是1G基因组的500k的时候 ,那一条contig 或者scaffold 的长度就叫做Contig N50和Scaffold N50.很明显这个数值越大说明组装的质量越好.

BACs, Fosmids。保持很好的共线性，说明拼接的很好。

Assembly Validation by cDNA dataset.

使用经典软件的优势是不用过多解释方法： broadinstitute

都用新软件则风险很大。方法review可能需要花很多时间。

最近institute发布了新分析流程。??//todo

end14:07

13. 重复序列去除

1.数据库repeatMasker(reference)只能识别已知重复； 2.repeatModeler(de novo)能识别未知重复。

end14:12

14. 基因组ncRNA

常规ncRNA、lncRNA、lincRNA、circRNA、snoRNA、piRNA
- Rfam: http://xfam.org/
miRNA: http://www.mirbase.org/
rRNA

15. 基因组组装

怎么挖掘 genome SSR marker ？
- 方便育种。方便被引用。
大片段扩增和基因重复：Lastz、Nucmer软件，进化分析分析。
trinity不好，但是无参拼接没有更好的了。所以尽量使用长reads。公司拼接可能去掉200bp以下的短reads，这样碎片少一点。
half de novo 就是使用近缘物种genome作为refer，拼接新物种。
protein map to genome: Spaln
PASA 真核生物基因组注释工具，用于refine。

end14:30

16. 基因开放阅读框/基因结构分析识别工具。

(部分基于预测的？因为有些不表达)

17. 基因组注释

没有近缘株，只用blast nr注释肯定不行。

需要做domain、motif分析。[P11,3]

hypothetical基本不算注释。
亮氨酸拉链 => 差不多是转录因子TF。

基因组比较：

寄生虫可能丢失脂肪酸代谢基因。

基因组也必须做KEGG分析：

原因：识别物种哪些代谢途径是完整的。
转录组表示在用该通路，基因组表示是否有该通路。
KEGG母图是不具有物种差异的。识别这个物种的代谢途径。

软件：多是一代测序时代的产物。

end14:44

18. 转录组分析[P12,3]

1.https://github.com/trinityrnaseq/RNASeq_Trinity_Tuxedo_Workshop/wiki

2.Tuxedo Genome Guided Transcriptome Assembly Workshop

3.RNA测序2个方向：真核转录组、RNA病毒测序。

4.一代、芯片关注：有polyA尾巴的RNA。相当数量lncRNA也有polyA，不编码蛋白。

二代可以关注repeat、假基因等。

ENCODE 计划：针对人和小鼠，完善基因组注释。https://www.encodeproject.org/data/annotations/

2008-09芯片、NGS较量厉害的2年。

需要3个生物学重复，不需要技术重复。

不要做等分取样？？？//todo

如果分析低丰度转录本，则加大数据量。
mRNA内标的问题。
- 至今没有被广泛接受的内标。
SAGE:1994-1995出来的技术。
- 开始以为能一次检测多个基因表达量，很受期待。
- 基因表达连续分析--SAGE(Serial Analysis of Gene Expression)
- 很难做，反复地PCR、电泳、回收，且几十bp片段很难回收。很鸡肋。
- refer http://www.dxy.cn/bbs/topic/266695?ppg=1
- 改进：增加片段长度。
- 实验很难，分析很难。
- small RNA测序时代：分离SAGE后，用小RNA盒子测序。(hASHIMOTO s-I,ET AL.(2009)PLoS ONE)
rRNA删除：
- 随着nonCoding的重视，开始从全RNA中删除rRNA。
- 删18S、28S序列：探针池，互相overlap覆盖18S+28S全长的一系列探针。增加捕获的效率。
- 做两轮删除，从80%下降到0.2%。
- rRNA删除试剂盒：人、拟南芥、G+、G-菌。
- 如果没有kit，先调取18S和28S，和人、拟南芥比对。自己合成一个盒子太贵！或者用人的捕获一次，用动物的捕获一次。
- 古菌：差别太大，很难做统一适用的盒子。clustalW看差别很大。简单粗暴：提取total RNA，直接测序。
适用性
- RNA | 研究明确，性状明确 |
- SAGE | eGWAS数量性状(红花+白花=>红、白、红花白点、白花红点，做co-expression，找表达差异基因模块，叫做eQTL) | 需要很多RNAseq。
- rRNA | 复杂调控机理 | 调控网络很多涉及很多ncRNA。|优点：分析全面。适合模式生物。|缺点：要求有高组装质量的基因组。成本最高的，>20G数据量。生信分析的压力最大。
lincRNA 功能很难验证。不编码蛋白，选择压力小。水稻xx基因，小麦的xx基因再做已经没有新意了，但是做lincRNA就很有新意，但实验难度并没有增加。
镇住销售的名词：核糖体去除链特异的RNAseq。
web分析系统：http://waprna.big.ac.cn/

19. RNAseq的常见图

5'3'完整性。
FPKM比RPKM分辨率高。[P16,1]
关心：中高表达、高表达的基因，基本都被研究烂了，是做质控用的，在文献中占用2-3行。
heatmap热图：可以做全基因的，DE基因，pathway基因等，单向、双向聚。分析发掘很重要。
- 小鼠1w、4w、10w的脑基因表达无基因symbol的热图分析：一周和十周基本呈现开和关状态，青春期小鼠实现从幼年到老年的关键阶段。
- 1w独有基因拿出来做GO分析，说明1w重点发育xx；
- 10w独有基因拿出来做GO分析，说明10w中xx在做什么。
芯片和表达的关系。中高表达量上一致，低表达量上NGS分辨率更高。芯片去除背景后，损失很多低表达基因。
生信需要学做图。基本是用R语言。http://www.genome.jp/kegg/ 中根据表达量显示颜色。
高表达、中低表达分析套路不同。
- 高丰度表达基因，DESeq等基于p值的，高表达一般是下游，集中在具体通路。
- 上游调控TF等低丰度基因，差异很小，p值也不显著。不用看P值，只用RPKM值相除，处理前/处理后，上调就染红，下调就染绿，给出整体的KEGG图也很有说服力。
- DAVID 做功能富集分析。橡胶割一刀后，分析到有损伤修复。
- IPA构建调控网络的构建（）B细胞分化过程相关基因调控网络。解读方式和DAVID不一样，一个球一个功能。网络的凝聚性，判断效应个数、关联。
直系基因（不同物种）、旁系基因（同一个物种扩增出来的基因）

20170326 day3

1.数据分析哟啊解决的问题？

mapping
转录组重建
表达水平定量
差异表达
图形化
- CummeRbund

2.fastqc质控

[P48]

质控后数据有问题怎么办？

lab5a 使用fastX Toolkit对数据进行清洗。[P74]

Available Tools

Fastq_to_fasta -h

FASTQ-to-FASTA converter
- Convert FASTQ files to FASTA files.
FASTQ Information
- Chart Quality Statistics and Nucleotide Distribution
FASTQ/A Collapser 去掉adaptor？
- Collapsing identical sequences in a FASTQ/A file into a single sequence (while maintaining reads counts)
FASTQ/A Trimmer
- Shortening reads in a FASTQ or FASTQ files (removing barcodes or noise).
FASTQ/A Renamer
- Renames the sequence identifiers in FASTQ/A file.
FASTQ/A Clipper
- Removing sequencing adapters / linkers
FASTQ/A Reverse-Complement
- Producing the Reverse-complement of each sequence in a FASTQ/FASTA file.
FASTQ/A Barcode splitter
- Splitting a FASTQ/FASTA files containning multiple samples
FASTA Formatter
- changes the width of sequences line in a FASTA file
FASTA Nucleotide Changer
- Convets FASTA sequences from/to RNA/DNA
FASTQ Quality Filter
- Filters sequences based on quality
FASTQ Quality Trimmer
- Trims (cuts) sequences based on quality
FASTQ Masker
- Masks nucleotides with 'N' (or other character) based on quality

  user120@ipags002:~/my_projects/lab5a_fastxtoolkit$ fastq_quality_filter -h
  usage: fastq_quality_filter [-h] [-v] [-q N] [-p N] [-z] [-i INFILE] [-o OUTFILE]
  Part of FASTX Toolkit 0.0.13 by A. Gordon (gordon@cshl.edu)

     [-h]         = This helpful help screen.
     [-q N]       = Minimum quality score to keep.
     [-p N]       = Minimum percent of bases that must have [-q] quality.
     [-z]         = Compress output with GZIP.
     [-i INFILE]  = FASTA/Q input file. default is STDIN.
     [-o OUTFILE] = FASTA/Q output file. default is STDOUT.
     [-v]         = Verbose - report number of sequences.
                    If [-o] is specified,  report will be printed to STDOUT.
                    If [-o] is not specified (and output goes to STDOUT),
                    report will be printed to STDERR.

  user120@ipags002:~/my_projects/lab5a_fastxtoolkit$

  #数据过滤
  $ fastq_quality_filter -q 30 -p 100 -i test_1.fq -o test_1_filter.fq -Q 33
  $ fastq_quality_filter -q 30 -p 100 -i test_2.fq -o test_2_filter.fq -Q 33

  #去掉前10个碱基
  $ fastx_trimmer -f 11 -i test_1_filter.fq -o test_1_quality_filter_trimmed.fq
  $ fastx_trimmer -f 11 -i test_2_filter.fq -o test_2_quality_filter_trimmed.fq

  #将序列小于50的序列去除
  $ fastx_clipper -l 50 -i test_1_quality_filter_trimmed.fq  -o test_1_quality_filter_trimmed_clipper.fq
  $ fastx_clipper -l 50 -i test_2_quality_filter_trimmed.fq  -o test_2_quality_filter_trimmed_clipper.fq

  # 数据格式转换
  #fastq2fasta格式
  $ fastq_to_fasta -i test_1.fq -o test_1.fasta

  #fasta2tabulated格式
  $ fasta_formatter -t -i test_1.fasta -o test_1_tabulated.txt

$ fastqc test_1_quality_filter_trimmed_clipper.fq

下载到本地,用浏览器查看fastqc结果，看数据清洗后数据质量是否提高了。

map method

un-spliced method:mapping to cDNA for quantification
- xx
spliced method: for gene structure construction
- Exon-first:适合绝大多数都能map。
- 局限：如果有含有intron的假基因，则可能直接map到假基因上了。
- seed and extend:打碎reads成小片段。
- 适用于unmap部分，或者unmap比例较高的转录组。

SAM文件格式

记录reads比对结果的文件格式。
SAM specs
- http://genome.sph.umich.edu/wiki/SAM
- https://en.wikipedia.org/wiki/SAM_(file_format)
- 11 列的纯文本tab文件。
- @开头是描述信息：specs文件后面有详述。
- SAM压缩成二进制后是BAM格式。

Col Field Type Regexp/Range Brief description
1 QNAME String [!-?A-~]{1,254} Query template NAME
2 FLAG Int [0,216-1] bitwise FLAG
3 RNAME String \*|[!-()+-<>-~][!-~]* Reference sequence NAME
4 POS Int [0,231-1] 1-based leftmost mapping POSition
5 MAPQ Int [0,28
-1] MAPping Quality
6 CIGAR String \*|([0-9]+[MIDNSHPX=])+ CIGAR string
7 RNEXT String \*|=|[!-()+-<>-~][!-~]* Ref. name of the mate/next read
8 PNEXT Int [0,231-1] Position of the mate/next read
9 TLEN Int [-231+1,231-1] observed Template LENgth
10 SEQ String \*|[A-Za-z=.]+ segment SEQuence
11 QUAL String [!-~]+ ASCII of Phred-scaled base QUALity+33

实验5b内容：使用bowtie2+sam tools +Bamview进行mapping、修改和查看。[]

建立基因组的index
- bowtie2-build xx_ref
align single-end reads
- bowtie2 -x xx_ref -U
align paired-end reads
align reads with local alighnment option
转成bam文件
排序sort
对bam做index
使用bamview查看结果。

lab5b：map到reference上

1.建立索引
$bowtie2-build ../../../datamart/sample_data/lab5b_ref.fa lab5b_ref
2.单端比对
$ bowtie2 -x lab5b_ref -U ../../../datamart/sample_data/lab5b_reads_1.fq -S eg1.sam
10000 reads; of these:
  10000 (100.00%) were unpaired; of these:
    596 (5.96%) aligned 0 times
    9404 (94.04%) aligned exactly 1 time
    0 (0.00%) aligned >1 times
94.04% overall alignment rate
user120@ipags002:~/my_projects/lab5b_bowtie2$

3.双端比对
$ bowtie2 -x lab5b_ref -1 ../../../datamart/sample_data/lab5b_reads_1.fq -2 ../../../datamart/sample_data/lab5b_reads_2.fq -S eg2.sam
10000 reads; of these:
  10000 (100.00%) were paired; of these:
    834 (8.34%) aligned concordantly 0 times
    9166 (91.66%) aligned concordantly exactly 1 time
    0 (0.00%) aligned concordantly >1 times
    ----
    834 pairs aligned concordantly 0 times; of these:
      42 (5.04%) aligned discordantly 1 time
    ----
    792 pairs aligned 0 times concordantly or discordantly; of these:
      1584 mates make up the pairs; of these:
        1005 (63.45%) aligned 0 times
        579 (36.55%) aligned exactly 1 time
        0 (0.00%) aligned >1 times
94.97% overall alignment rate

4.使用本地比对选项。
$ bowtie2 --local -x lab5b_ref -U ../../../datamart/sample_data/longreads.fq -S eg3.sam
6000 reads; of these:
  6000 (100.00%) were unpaired; of these:
    157 (2.62%) aligned 0 times
    5634 (93.90%) aligned exactly 1 time
    209 (3.48%) aligned >1 times
97.38% overall alignment rate

5.sam2bam
$ samtools view -bS eg2.sam > eg2.bam

6.bam排序
samtools sort eg2.bam eg2.sorted

7.为排序过的bam建lindex
samtools index eg2.sorted.bam

8.查看比对结果
mkdir eg2
cp ../../../datamart/sample_data/lab5b_ref.fa eg2.ref.fa
cp eg2.ref.fa eg2.sorted.bam eg2.sorted.bam.bai eg2
下载eg2文加夹。
打开bamview文件，输入bam文件和eg2.ref.fa查看比对结果。

Bamview接收输入bam文件，可选输入reference。

蓝色同向
黑色反向
红色
可以一直放大，直到看到碱基。
可以单击序列，查看配对序列位置。

更好的是使用IGV可视化。

转录组重建算法

De brijn graph
- k-mer is a substring of length K.
- 3mer, 分成2mer，left和right。
- 然后序列按照mer进行画线连接。
- 解析这些图的算法的实现：
- scirpture: 选择所有路径？？
- cufflink: 选择够用的路径？？

基因组依赖的组装

两个补充
- 分清楚错误和变异
- 灵敏度和复杂度最优平衡。
去掉low-covverage分支。
k-mer length->不同结果。
- AbYss -> variable k-mer strategy.

refer links: http://www.htslib.org/doc/

2017.3.26

1. 基因表达定量

FPKM(用于双端PE测序) vs RPKM
转录水平
- Unique reads only method 低估一些低表达isoform
- isoform level 定量 using MLE(most likely estimateds)[best]
gene level
- map到exon的reads 计数
- map到
- ？
有参：
- 参考基因组
- 参考注释文件
  - IGV可以加track，文件格式是GFF格式

GFF3文件格式：9列tab文件，记录基因信息。

Column 1: "seqid"
Column 2: "source"
Column 3: "type"
Columns 4 & 5: "start" and "end"
Column 6: "score"
Column 7: "strand"
Column 8: "phase"
Column 9: "attributes" 预定义很多标签(id,parent,)

0  ##gff-version 3.2.1
 1  ##sequence-region ctg123 1 1497228
 2  ctg123 . gene            1000  9000  .  +  .  ID=gene00001;Name=EDEN
 3  ctg123 . TF_binding_site 1000  1012  .  +  .  ID=tfbs00001;Parent=gene00001
 4  ctg123 . mRNA            1050  9000  .  +  .  ID=mRNA00001;Parent=gene00001;Name=EDEN.1
 5  ctg123 . mRNA            1050  9000  .  +  .  ID=mRNA00002;Parent=gene00001;Name=EDEN.2
 6  ctg123 . mRNA            1300  9000  .  +  .  ID=mRNA00003;Parent=gene00001;Name=EDEN.3
 7  ctg123 . exon            1300  1500  .  +  .  ID=exon00001;Parent=mRNA00003
 8  ctg123 . exon            1050  1500  .  +  .  ID=exon00002;Parent=mRNA00001,mRNA00002
 9  ctg123 . exon            3000  3902  .  +  .  ID=exon00003;Parent=mRNA00001,mRNA00003
10  ctg123 . exon            5000  5500  .  +  .  ID=exon00004;Parent=mRNA00001,mRNA00002,mRNA00003
11  ctg123 . exon            7000  9000  .  +  .  ID=exon00005;Parent=mRNA00001,mRNA00002,mRNA00003
12  ctg123 . CDS             1201  1500  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1
13  ctg123 . CDS             3000  3902  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1

能从gff文件看出复杂的基因结构的注释信息：

一个基因id有3个mRNA
有4个外显子
可以导入gff文件到IGV，生成track。
其他基因注释格式
基因组文件 + gff文件 => 导出转录本信息、promoter信息等。
更多基因相关文件格式：https://genome.ucsc.edu/FAQ/FAQformat.html
bed文件时gff的简化版本。一般都是tab文件。

2.差异表达分析算法

DE 分析方法。
不同condition的mean sv 画分布图。
modeling the distribution of count.
- 泊松分布可能最合适，后来又发展其他non-linear方法。
基因定量的3种方法。
- 直接map到exon
- splice map 到exon
- isoform更好
考虑的其他因素
- 测序深度
- 基因表达量
- 基因长度
- lib构建方法
- 只有生物学重复才能提高严谨度。
数据分析软件
- 有参考*
- 有基因组序列
- 有基因组注释信息
- Tuxedo(燕尾服) package
lab6 有参RNAseq分析
1. bowtie、 tophat
2. cufflink
3. cummeRbund
4. 如果不想发现新的选择性剪切，tophat加参数 -no-??
5. 推荐：每个样品分别做cufflink，在cuffmerge合并。
6. 要有对照，要有重复。怎么描述？TSS I有2个转录本，TSS II有1个转录本。哪个基因转录因子用的多？哪个转录本多？

3.实验

$ ls 4k_READS_sample cummeRbund.R.script genes.gtf genome.fa
http://agetouch.blog.163.com/blog/static/2285350902016739843185/

1.建立索引
$ bowtie2-build genome.fa genome

2.tophat 做mapping
$ tophat -p 4 -G genes.gtf -o C1_R1_thout genome 4k_READS_sample/C1_R1_1.fq 4k_READS_sample/C1_R1_2.fq
$ tophat -p 4 -G genes.gtf -o C1_R2_thout genome 4k_READS_sample/C1_R2_1.fq 4k_READS_sample/C1_R2_2.fq

$ tophat -p 4 -G genes.gtf -o C2_R1_thout genome 4k_READS_sample/C2_R1_1.fq 4k_READS_sample/C2_R1_2.fq
$ tophat -p 4 -G genes.gtf -o C2_R2_thout genome 4k_READS_sample/C2_R2_1.fq 4k_READS_sample/C2_R2_2.fq

生成bam文件。最重要的是accepted_hits.bam文件，在cufflinks和cuffdiff都会用到。

3.组装基因和转录组
$ cufflinks -p 4 -o C1_R1_clout C1_R1_thout/accepted_hits.bam
$ cufflinks -p 4 -o C1_R2_clout C1_R2_thout/accepted_hits.bam
$ cufflinks -p 4 -o C2_R1_clout C2_R1_thout/accepted_hits.bam
$ cufflinks -p 4 -o C2_R2_clout C2_R2_thout/accepted_hits.bam

生成gtf文件，transcripts.gtf。
会有基因定量，gene.tracking,isoxx.tracking文件。

echo "./C1_R1_clout/transcripts.gtf" > assemblies.txt
echo "./C1_R2_clout/transcripts.gtf" >> assemblies.txt
echo "./C2_R1_clout/transcripts.gtf" >> assemblies.txt
echo "./C2_R2_clout/transcripts.gtf" >> assemblies.txt

4.创建单个的转录注释文件
cuffmerge -g genes.gtf -s genome.fa -p 4 assemblies.txt

5. 识别差异表达基因和转录本
cuffdiff -o diff_out -b genome.fa -p 4 -L C1,C2 -u merged_asm/merged.gtf ./C1_R1_thout/accepted_hits.bam,./C1_R2_thout/accepted_hits.bam ./C2_R1_thout/accepted_hits.bam,./C2_R2_thout/accepted_hits.bam

注意2个条件之间要有空格，条件内部是没有空格的逗号。

exp.diff是最重要的文件。

MA-Plot怎么看？

http://bioinfo.cipf.es/babelomicstutorial/maplot

MA图用于研究2个变量的log比值，和这两个变量的的均值的依赖性。
Y轴是M值（log尺度的minus）是2个变量的log比值。
X轴是A值，2个参数的平均值（log尺度的average）。
X: logCounts
Y: logFc.
应该是一条线？？

volcano plot

FDR
- 考虑多样本矫正，又考虑?的排名。
- 做法：p-value * 测序基因数量 / ?排名。

无参转录组

非模式生物数据分析软件。

trinity(三位一体) software

Inchworm: 构建contigs，大概是isoform
Chrysalis cluster: contig clusters. de Bruijn graphs
Butterfly: reconstructed isoforms

trinity程序的组成部分
- reads组装需要大型服务器
- 聚类：连接contig成cluster
- de brujn graph，之后解析基因等。

运行trinity的硬件：512G-1T内存。
trinity: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3571712/
https://github.com/trinityrnaseq/trinityrnaseq/wiki

使用UCSC基因组浏览器看ENCODE数据：https://www.genome.gov/multimedia/slides/encode2015-researchappsusers/08_fujita.pdf

20170326 Linux命令

$ ssh user120@123.57.214.96
961251

视频链接在ppt中。 QQ群: 554135121

概念与历史

2000 fly

WGS
生物学意义 biology of genome
疾病 biology of disease
制药 medicine
个性化医疗 healthcare
Nextera:打断DNA。
为什么Illunima测序不是很长？
PE测序顺序：read1,index1,reverse,index2,PCR, read2.
mate pair 和pair end的区别？
- matePair建库不同。两头变中间，中间变两头。相当于中间测得更长了。
三代测序比较成熟的pacBio和nanoPore
pacBio的single molecule seq的2项技术:
1. 核酸磷酸化连接（以前是荧光连接到碱基上，现在是连接到磷酸上。）
2. 纳米小孔底部固定DNA聚合酶，高信噪比光学探头聚焦到孔底部。还要快速成像，很小的孔。荧光加上和解开时间差别大。
nanoPore的
1. 正反链都测
2. 使用一个膜蛋白，放到一个不导电膜上，监测电流，如果有分子走过去，电流就会降低。

中草药参考基因库（植物所）。

RNAseq实验设计

什么能发高分？

discovery 新机理
resource 测序了新的基因组
method 发明了新的方法
tools 网站、软件

RNA-seq的目的：

表达量
差异表达

RNA-seq和microarray的比较

中表达基本一致
低表达、高表达不一致。

测序深度的选择

深度越高，则发现的基因越多。
但是多高到一定程度，就不会发现新基因了。
fasta格式.
quality score: Q=-log10(P). Q+33的ASCII码描述碱基，得到fastq格式。

fastqc数据质控

1.文件权限

chown UID:GID files chmod o+w file.txt #other组增加write权限

文本处理

grep split swk

文件查找

find

sh script

! /bin/bash

date ls -la

文本编辑

nano emacs vi

进程管理

su suer passwd top top -u user history ps ps -ef 显示所有账号的进程详细列表 kill -9 pid

几个特殊符号

< 导入符号到文件 | 管道符号，输出作为下一步的输入

重定向 & 后台处理

软件安装

ubuntu 使用apt-get red hat使用 yum

加压缩

tar -xzvf xx.tar.gz ./configure make (sudo) make install make clean

注意

1.注意大小写

end8:34

fastqc

http://yanshouyu.blog.163.com/ http://yanshouyu.blog.163.com/blog/static/214283182201302835744453/

http://cbsu.tc.cornell.edu/lab/doc/RNASeq_workshop_2013_part2.pdf

2017.3.27

基因功能注释

通过表达量FC，p-value，得到DE gene list。
GO & KEGG。
trinotate 转录本功能注释
- blast
- pfam 蛋白注释
- uniprot
- SQLite

trinotate流程

TransDecoder 预测编码区。
blast获取同源序列。
HMMER 使用隐马尔科夫模型发现蛋白质domain，很敏感。
signalP 发现信号肽。
tmHMM 预测跨膜区。
RNAMMER 识别rRNA。
Trinotate 会输入如上结果到SQLite数据库。
以上1-6的server服务需要自己运行，之后把结果放到一个文件夹下，给步骤7，输出为excel文件。

转录组基因注释
lab8a

lab8b

差异基因列表的目的

列表中什么基因？
list中的基因家族是什么？
富集的term是什么？
？？

什么是好的term？

包含很多重要的基因（marker gene）
合适的基因数（100-2000）
大多数基因统计上显著
上下调基因有生物学意义
有富集生物学意义
高度重复
其他wet实验可以验证的。

CASEG 分析

Cluster of Adjacent and Similarly Expressed Genes(CASEG)

https://www.researchgate.net/publication/221586673_The_Human_Genome_Contains_Numerous_Clusters_of_Adjacent_Co-Regulated_Genes

GO一级分类3个
数据预处理 log之后可能就正态分布了。
富集分析也是很复杂的领域，剩下的不过：DAVID等。
- class I: singular enrichment Analysis仅仅考虑p值。
- class II:考虑基因的FC和p-value。GSEA。
- class III:还要考虑DE gene之间的关系。wgcna。
- Gene命名产生了GO，后来GO越来越复杂了。Gene Oncology.
- GO enrichment analysis
- map gene 到Goterm ，blast2GO
- GO 画图。
- 蛋白互作网络PPI: string-db,
- gene list analysis: GoSeq, KOBAS, DAVID[推荐].
- 非模式生物，不建议做GO分析，没有意义。
- DAVID最早使用affimatrix的probe ID作为输入的，现在支持的输入越来越多了。
http://geneontology.org/page/go-enrichment-analysis
- AltAnalyze
- for model organism
- was designed for ??

全长转录本测序技术IsoSeq

模式生物等已经研究的差不多了，就不用做三代测序了。三代不能定量。仅用于发现新长转录本。
前提是能够用实验手段获得全长转录本。
样本制备时还是有PCR。目前测序过程中酶能跑到30k。
建库：bluePipin仪器或者手工切胶建库。
Hybrid assembly is a frontier.

DawnEve / bioToolKit

2017-"高通量测序应用最新技术与数据分析"培训班笔记 #1