eigen 列拼接_cufflinks

短命女 2023-01-01 06:48 192阅读 0赞

1. 数据

手头有两个转录组数据，分别是某一真菌物种(该物种没有基因组数据)的双核菌丝阶段和子实体阶段的转录组数据。测序平台是Illumina Hiseq2000；插入片段长度200bp，测序的reads长度90bp。得到的数据文件为：

/home/user/RNA-seq/mycelium\_reads1.fastq

/home/user/RNA-seq/mycelium\_reads2.fastq

/home/user/RNA-seq/fruitingbody\_reads1.fastq

/home/user/RNA-seq/fruitingbody\_reads2.fastq

2. 数据的预处理

2.1. 去除N的比例大于5%的reads;去除低质量reads(质量值Q≤20的碱基数占整个read的50％以上)；

2.2. 根据FastQC对上述过滤后的reads的质量检测，去除reads首尾各10bp的碱基，得到的预处理数据为：

/home/user/RNA-seq/clean\_reads/mycelium\_reads1.fastq

/home/user/RNA-seq/clean\_reads/mycelium\_reads2.fastq

/home/user/RNA-seq/clean\_reads/fruitingbody\_reads1.fastq

/home/user/RNA-seq/clean\_reads/fruitingbody\_reads2.fastq

3. 使用Trinity和TGICL进行转录组的组装

3.1. 对mycelium数据和fruitingbody数据分别进行转录组的组装

$ pwd

/home/user/RNA-seq/

$ mkdir assembly

$ cd assembly

$ Trinity.pl --jaccard\_clip --seqType fq --JM 50G --SS\_lib\_type FR --CPU 24 --inchworm\_cpu 24 --bflyCPU 24 --group\_pairs\_distance 500 -min\_contig\_length 200 --output mycelium\_contig --left /home/user/RNA-seq/clean\_reads/mycelium\_reads1.fastq --right /home/user/RNA-seq/clean\_reads/mycelium\_reads2.fastq

$ Trinity.pl --jaccard\_clip --seqType fq --JM 50G --SS\_lib\_type FR --CPU 24 --inchworm\_cpu 24 --bflyCPU 24 --group\_pairs\_distance 500 -min\_contig\_length 200 --output fruitingbody\_contig --left /home/user/RNA-seq/clean\_reads/fruitingbody\_reads1.fastq --right /home/user/RNA-seq/clean\_reads/fruitingbody\_reads2.fastq

当然，需要将生成的两个转录组的Trinity.fasta序列按长度进行排序；统一更改的fasta

头名称；更改fasta文件名

$ cp mycelium\_contig/Trinity.fasta mycelium\_contigs.fasta

$ cp fruitingbody\_contig/Trinity.fasta fruitingbody\_contigs.fasta

3.2. 使用TGICL将两个转录组序列合并

$ pwd

/home/user/RNA-seq/assembly

$ mkdir all\_tissue\_contig

$ cd all\_tissue\_contig

$ cat ../mycelium\_contigs.fasta ../fruitingbody\_contigs.fasta > all.contigs.fasta

$ tgicl -F all.contigs.fasta

tgicl生成了两个有用的文件： asm\_1/contigs 和 all.contigs.fasta.single

tons。其中前者是聚类后的contigs结果；后者是没有聚类的单独的contigs的序列名，需

要分别到../mycelium\_contigs.fasta 和 ../fruitingbody\_contigs.fasta文

件中提取出相应的序列: all.contigs.fasta.singletons.mycelium.fasta 和

all.contigs.fasta.singletons.fruitingbody.fasta

$ cat asm\_1/contigs all.contigs.fasta.singletons.mycelium.fasta all.contigs.fasta.singletons.fruitingbody.fasta > all\_contigs.fasta

在对all\_contigs.fasta进行fasta头的重命名，并将序列按长度排序；同时要得到all\_

contigs.fasta中的序列和mycelium\_contigs.fasta，fruitingbody\_contigs.

fasta中序列的对应关系。

$ cp all\_contigs.fasta ../

3.3. 至此得出转录组的组装结果：

/home/user/RNA-seq/assembly/all\_contigs.fasta

4. 使用cufflinks来分析差异表达基因

4.1 使用tophat来将预处理后的reads比对到转录组序列上

$ pwd

/home/user/RNA-seq/

$ mkdir cufflinks

$ cd cufflinks

$ bowtie2-build ../all\_contigs.fasta all\_contigs

$ tophat -o tophat\_out\_mycelium -r 60 --mate-std-dev 80 -p 24 --library-type fr-unstranded all.contigs /home/user/RNA-seq/clean\_reads/mycelium\_reads1.fastq /home/user/RNA-seq/clean\_reads/mycelium\_reads2.fastq

$ tophat -o tophat\_out\_fruitingbody -r 60 --mate-std-dev 80 -p 24 --library-type fr-unstranded all.contigs /home/user/RNA-seq/clean\_reads/fruitingbody\_reads1.fastq /home/user/RNA-seq/clean\_reads/fruitingbody\_reads2.fastq

4.2 自制一个transcripts.gtf文件

由于cuffdiff运行需要一个参考序列的transcripts.gtf文件: 该文件有9列，使用tab分隔；使exon的范围为整个contig。其格式如：

AA.auricula\_all\_contig\_1 chenlianfu exon 1 9401 . . . gene\_id "A.auricula\_all\_contig\_1"; transcript\_id "A.auricula\_all\_contig\_1";

A.auricula\_all\_contig\_10 chenlianfu exon 1 4464 . . . gene\_id "A.auricula\_all\_contig\_10"; transcript\_id "A.auricula\_all\_contig\_10";

A.auricula\_all\_contig\_100 chenlianfu exon 1 3090 . . . gene\_id "A.auricula\_all\_contig\_100"; transcript\_id "A.auricula\_all\_contig\_100";

A.auricula\_all\_contig\_1000 chenlianfu exon 1 1768 . . . gene\_id "A.auricula\_all\_contig\_1000"; transcript\_id "A.auricula\_all\_contig\_1000";

A.auricula\_all\_contig\_10000 chenlianfu exon 1 586 . . . gene\_id "A.auricula\_all\_contig\_10000"; transcript\_id "A.auricula\_all\_contig\_10000";

A.auricula\_all\_contig\_10001 chenlianfu exon 1 586 . . . gene\_id "A.auricula\_all\_contig\_10001"; transcript\_id "A.auricula\_all\_contig\_10001";

4.3 使用cuffdiff来分析转录子的表达量和差异表达基因

$ pwd

/home/user/RNA-seq/cufflinks

$ cuffdiff -L mycelium,fruitingbody --library-type fr-unstranded -p 8 -o cuffdiff ./transcriptome.gtf ./tophat\_out\_mycelium/accepted\_hits.bam ./tophat\_out\_fruitingbody/accepted\_hits.bam

至此得出转录子的表达量数据和差异表达分析

/home/user/RNA-seq/cufflinks/cuffdiff/gene\_exp.diff