搜索
查看: 636|回复: 4

有参转录组分析流程 -- 仅供参考

[复制链接]

1

主题

1

帖子

24

积分

版主

Rank: 7Rank: 7Rank: 7

积分
24
发表于 2016-9-19 17:55:06 | 显示全部楼层 |阅读模式
    从0开始,先分享几个学习资源:
  • 生信入门视频等教程
    http://pan.baidu.com/s/1kUq1xLP 密码:i829
    http://pan.baidu.com/s/1ntUIBix 密码:lluw
  • 网站1(Bioinformatics Team (BioITeam) at the University of Texas)
    https://wikis.utexas.edu/display/bioiteam/Differential+gene+expression+analysis
  • 生信群(生信菜鸟团)
    http://www.bio-info-trainee.com/?p%3D771
  • 网站2
    http://bioinformatics.ca/
  • 博客(糗哥的文)
    http://blog.qiubio.com:8080/archives/3007/comment-page-2


      上述资源挑点看看应该能明白生信是干啥的了吧(我可是讲不明白的,汗),剩下的就是实战喽。还是要找数据过一遍才能发现问题,解决问题。出问题了真是一种折磨呀,但是,痛苦了才会有进步嘛。。。。。。
      下面是有参转录组分析流程:

步骤:
1. fastqc()
2. NGSQC()
3. bowtie(建立索引)
4. tophat(比对)
5.1 RSeQC (RNA-seq 数据质控)            
5.2 cufflinks  cuffmerge cuffdiff(拼接组装,寻找差异基因)  cummeRbund(可视化)
6. HTseq(基因表达量计算)
7. DESeq                  or           edgeR (寻找差异基因)
8. GO  (GO和KEGG富集的R包clusterProfiler)             or           KEGG
    GeneOntology                                                           or           kobas
10. TopGo(富集分析可视化)

一、原始数据测序质量分析
raw_data(大豆-Soybean,42G)

1.1 fastQC
输出文件:fastqc-report            
制作接头文件:
将Overrepresented sequences中的重复序列放在接头文件中,以>universal开头




1.2 NGSQC
命令行:
         perl /home/soft/out_softwares/NGSQCToolkit_v2.3.3/QC/IlluQC_PRLL.pl -pe ../reads/s3_ZH10_72h_R1.fastq ../reads/s3_ZH10_72h_R2.fastq adaptor_s3.txt A

         perl /home/soft/out_softwares/NGSQCToolkit_v2.3.3/QC/IlluQC_PRLL.pl -pe ../reads/s7_BF15_72h_R1.fastq ../reads/s7_BF15_72h_R2.fastq adaptor_s7.txt



二、 Mapping至参考基因组
2.1 bowtie
bowtie2-build ../ref/Gmax_275_v2.0.fa ../ref/Gmax_275_v2.0
2.2 tophat
tophat2 -o S3-1-thout ../ref/Gmax_275_v2.0 ../reads/IlluQC_Filtered_files/s3_ZH10_72h_R1.fastq_filtered ../reads/IlluQC_Filtered_files/s3_ZH10_72h_R2.fastq_filtered

tophat2 -o S7-1-thout ../ref/Gmax_275_v2.0 ../reads/IlluQC_Filtered_files/s7_BF15_72h_R1.fastq_filtered ../reads/IlluQC_Filtered_files/s7_BF15_72h_R2.fastq_filtered



三、数据质控
         首先查看bam文件是否sort:  samtools view -h ../tophat/S3-1-thout/accepted_hits.bam |head,如果出现SO:coordinate,则表示排序过。
3.1 RSeQC
bam_stat.py -i ../tophat/S3-1-thout/accepted_hits.bam>S3-1-bam_stat.txt
bam_stat.py -i ../tophat/S7-1-thout/accepted_hits.bam>S7-1-bam_stat.txt
clipping_profile.py -i ../tophat/S3-1-thout/accepted_hits.bam -o S3-1-clipping_profile-output
clipping_profile.py -i ../tophat/S7-1-thout/accepted_hits.bam -o S7-1-clipping_profile-output
mismatch_profile.py -i ../tophat/S3-1-thout/accepted_hits.bam -o S3-1-clipping_profile-output
mismatch_profile.py -i ../tophat/S7-1-thout/accepted_hits.bam -o S7-1-clipping_profile-output
insertion_profile.py -i ../tophat/S3-1-thout/accepted_hits.bam -o S3-1-insertion_profile-output
insertion_profile.py -i ../tophat/S7-1-thout/accepted_hits.bam -o S7-1-insertion_profile-output
bam2wig.py -i ../tophat/S3-1-thout/accepted_hits.bam -o S3-1-bam2wig-output
bam2fq.py -i ../tophat/S7-1-thout/accepted_hits.bam -o S7-1-bam2wig-output


四、表达差异分析
(一)

4.1 cufflinks
cufflinks -o S3-1-clout ../tophat/S3-1-thout/accepted_hits.bam
cufflinks -o S7-1-clout ../tophat/S7-1-thout/accepted_hits.bam

#touch assemblies.txt
echo "./E2-1-clout/transcripts.gtf" >> assemblies.txt
echo "./E2-2-clout/transcripts.gtf" >> assemblies.txt

cuffmerge -s ../ref/Gmax_275_v2.0.fa assemblies.txt
## 需要注释文件
cuffdiff -o diff_out1 -b ../ref/Gmax_275_v2.0.fa -L S3-1,S7-1 -u merged_asm/merged.gtf ../tophat/S3-1-thout/accepted_hits.bam ../tophat/S7-1-thout/accepted_hits.bam
cuffdiff -o diff_out2 -b ../ref/Gmax_275_v2.0.fa -L S3-1,S7-1 -u ../ref/Gmax_275_Wm82.a2.v1.gene.gff3 ../tophat/S3-1-thout/accepted_hits.bam ../tophat/S7-1-thout/accepted_hits.bam

4.2 cummeRbund包可视化
/data/mariel_ma/cummeRbund
/data/mariel_ma/cummeRbund/plot/cuff_plot.R
https://yunpan.cn/cxtgtkYk5iRak  访问密码 ee0f



(二)

4.1 HTseq
samtools view -h ../tophat/S3-1-thout/accepted_hits.bam >S3-1-accepted_hits.sam
samtools view -h ../tophat/S7-1-thout/accepted_hits.bam >S7-1-accepted_hits.sam

/opt/anaconda/bin/htseq-count -m intersection-nonempty -t gene -i ID ./S3-1-accepted_hits.sam ../ref/Gmax_275_Wm82.a2.v1.gene.gff3 >S3-1.sam.count

/opt/anaconda/bin/htseq-count -m intersection-nonempty -t gene -i ID ./S7-1-accepted_hits.sam ../ref/Gmax_275_Wm82.a2.v1.gene.gff3 >S7-1.sam.count

4.2 DEseq
join ../htseq/S3-1.sam.count ../htseq/S7-1.sam.count >gene_counts_HTseq.count

cat gene_counts_HTseq.count|sed 's/\s/\t/g'|sed '1i\locus_tag\tS3\tS7' > gene_counts.tab

perl /home/soft/out_softwares/Trinity/trinityrnaseq-2.0.6/Analysis/DifferentialExpression/run_DE_analysis.pl --matrix gene_counts.tab --method DESeq

4.3 edgeR

join ../htseq/S3-1.sam.count ../htseq/S7-1.sam.count >gene_counts_HTseq.count


cat gene_counts_HTseq.count|sed 's/\s/\t/g'|sed '1i\locus_tag\tS3\tS7' > gene_counts.tab

perl /home/soft/out_softwares/Trinity/trinityrnaseq-2.0.6/Analysis/DifferentialExpression/run_DE_analysis.pl --matrix gene_counts.tab --method edgeR




五、KEGG
(一)


5.1 kobas
cat gene_counts.tab.S3_vs_S7.edgeR.DE_results|sed '1d' |awk -F "\t" '{print $1}' >geneid.txt
perl trans_id.pl Gmax_275_geneid2entrez.txt  1 2 geneid.txt>kegg2id.txt
run_kobas.py -i kegg2id.txt -t id:ncbigene -s gmx -d K -o kegg_pathway.txt



本文数据: 数据太大了,传了几次没传上去,我再想想办法传上去。。。。。。


文中涉及文件:https://yunpan.cn/cxWcEV76zixJu  访问密码 4472

本文参考:
1. http://www.vcru.wisc.edu/simonlab/bioinformatics/programs/trinity/docs/analysis/diff_expression_analysis.html





上一篇:你的基因名称写对了吗
下一篇:基因组求助帖:platanus 高杂合基因组的组装
回复

使用道具 举报

0

主题

21

帖子

146

积分

注册会员

Rank: 2

积分
146
发表于 2016-11-2 21:22:38 | 显示全部楼层
对于有重复的样本,DEseq,edgeR怎么设计gene_counts.tab?
回复 支持 反对

使用道具 举报

0

主题

21

帖子

146

积分

注册会员

Rank: 2

积分
146
发表于 2016-11-2 22:18:59 | 显示全部楼层
云盘链接都无效
回复 支持 反对

使用道具 举报

2

主题

4

帖子

125

积分

注册会员

Rank: 2

积分
125
QQ
发表于 2016-12-27 11:39:52 | 显示全部楼层
云盘分享链接失效了,再来一次吧。
Yes, I am serious.
回复 支持 反对

使用道具 举报

1

主题

13

帖子

98

积分

注册会员

Rank: 2

积分
98
发表于 2017-1-10 15:04:12 | 显示全部楼层
raw data 数据处理部分用什么?
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|关于我们|手机版|小黑屋|生信技能树    

GMT+8, 2017-1-25 01:06 , Processed in 0.189553 second(s), 33 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.