搜索
查看: 11799|回复: 3

[mRNA-seq] 转录组入门(1-6)从测序数据到生成count矩阵

[复制链接]

29

主题

131

帖子

1208

积分

金牌会员

Rank: 6Rank: 6

积分
1208
发表于 2017-7-28 22:02:37 | 显示全部楼层 |阅读模式
本帖最后由 anlan 于 2017-7-28 22:15 编辑

文献名称:AKAP95 regulates splicing through scaffolding
RNAs and RNA processing factors
  • 查找数据:Data availability
    The RIP-seq an RNA-seq data have been deposited in the Gene
    Expression Omnibus database, with accession code GSE81916. All other data is
    available from the author upon reasonable request.
  • 获得GSE号:GSE81916
  • https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE81916获取数据信息,并点击网址下方的ftp,下载测序数据
  • https://trace.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA323422可知我们需要的mRNA测序编号为SRR3589956到SRR3589962
  • 通过Apera下载SRR数据,这里以SRR3589956为例:
    [AppleScript] 纯文本查看 复制代码
    ascp -T -i /home/anlan/.aspera/connect/etc/asperaweb_id_dsa.openssh [url=mailto:anonftp@ftp-private.ncbi.nlm.nih.gov]anonftp@ftp-private.ncbi.nlm.nih.gov[/url]:sra/sra-instant/reads/ByRun/sra/SRR/SRR358/SRR3589956/SRR3589956.sra ./
  • 通过sratoolkit工具将SRR文件转化为fastq格式的测序数据(写了个shell循环)
    [Shell] 纯文本查看 复制代码
    for i in $(seq 56 62);do nohup fastq-dump --split-3  SRR35899${i} &;done
    
  • 通过fastqc对每个fastq文件进行质检,用multiqc查看整体质检报告(对当前目录下的fastq测序结果进行质检,生成每个fq文件的质检报告总multiqc整合后统计查看)
    [AppleScript] 纯文本查看 复制代码
    fastqc *.fastq
    multiqc ./

    点击这个url可以查看我这个multiqc报告:http://www.bioinfo-scrounger.com/data/multiqc_report.html
  • 如果有接头或者质量值不达标的需要进行过滤,这次的数据质量都不错,因此直接进行比对即可
  • 安装hisat2软件,下载人类的hiast2索引文件

  • 用hisat2进行比对,测序数据放在data目录下,索引文件放在reference/index/hisat2/hg19目录下,SRR3589956-SRR3589958为人的测序数据
    [Shell] 纯文本查看 复制代码
    for i in $(seq 56 58);do hisat2 -p 4 \-x ~/reference/index/hisat2/hg19/genome \-1 ./data/SRR35899${i}_1.fastq -2 ./data/SRR35899${i}_2.fastq \-S SRR35899$i.sam >SRR35899${i}.log;done
  • 用samtools将sam文件转化为bam文件,并使用默认排序
    [Shell] 纯文本查看 复制代码
    for i in $(seq 56 58);do samtools sort -@ 5 -o SRR35899${i}.bam SRR35899${i}.sam;done
  • 用htseq对比对产生的bam进行count计数

    • htseq安装,使用miniconda,省事!唯一的问题是htseq版本不是最新的,是0.7.2。想要最新版还是要正常安装,可参考http://www.biotrainee.com/thread-1847-1-2.html
        
      [AppleScript] 纯文本查看 复制代码
      conda install -c bioconda htseq
    • 用htseq将对比后的结果进行计数
        
      [Shell] 纯文本查看 复制代码
      for i in $(seq 56 58);do htseq-count -f bam -r pos -s no \
                              SRR35899${i}.bam ~/reference/genome/hg19/gencode.v26lift37.annotation.gtf \
                              1>SRR35899${i}.count 2>SRR35899${i}_htseq.log;done

  • 将3个count文件(SRR3589956.count,SRR3589957.count,SRR3589958.count)合并成一个count矩阵,这是就需要脚本来解决这个问题,不然其他方法会稍微麻烦点
    [Shell] 纯文本查看 复制代码
    #!/usr/bin/perl -w
    use strict;
    my $path = shift @ARGV;
    opendir DIR, $path or die;
    my @dir = readdir DIR;
                    
    my $header;
    my @sample;
    my %hash;
    foreach my $file (@dir) {
            if ($file =~ /^\w+.*\.count/) {
                    push @sample, $file;
                    $header .= "\t$file";
                    open my $fh, $file or die;
                    while (<$fh>) {
                            chomp;
                            next if ($_ =~ /^\W+/);
                            my @array = split /\t/, $_;
                            $hash{$array[0]} -> {$file} = $array[1];
                    }
                    close $fh;
            }
    }
    print "$header\n";
    map{
            my $gene = $_;
            print "$gene";
            foreach my $file (@sample) {
                    print "\t".$hash{$gene} -> {$file};
            }
            print "\n";
    }keys %hash;
  • 按照接下来的剧本,应该讲count_matrix文件导入DESeq进行差异表达分析。但是从这篇文章的Bioinformatic analyses部分可以发现,作者的control组的2组数据是来自2个不同的批次(一个是SRR3589956,另外一个来源GSM1095127 in GSE44976),treat组倒是同一个批次(SRR3589957和SRR3589958)。但是对于Mouse cells来说,倒是满足2个control和2个treat都正常来自同个批次,因此打算重新用SRR3589959-SRR3589962重新做个一个count_matrix进行后续差异分析

有道笔记链接:http://note.youdao.com/noteshare ... 09b4f004d8f4542e36e





上一篇:【高薪招聘】生物信息科学家/工程师
下一篇:(0)系列如何入门生物信息学
回复

使用道具 举报

3

主题

25

帖子

464

积分

中级会员

Rank: 3Rank: 3

积分
464
发表于 2017-9-15 12:18:26 | 显示全部楼层
厉害厉害,佩服佩服!
回复 支持 反对

使用道具 举报

0

主题

5

帖子

71

积分

注册会员

Rank: 2

积分
71
发表于 2018-11-14 09:57:20 | 显示全部楼层
[root@dodogene etc]# ascp -T -i /root/.aspera/connect/etc/asperaweb_id_dsa.openssh anonftp@ftp-private.ncbi.nlm.nih.gov :sra/sra-instant/reads/ByRun/sra/SRR/SRR358/SRR3589956/SRR3589956.sra ./
ascp: no remote host specified, exiting.
回复 支持 反对

使用道具 举报

0

主题

5

帖子

71

积分

注册会员

Rank: 2

积分
71
发表于 2018-11-14 10:52:42 | 显示全部楼层
数据下载好了,不过转化时候,for i in $(seq 56 62);do nohup fastq-dump --split-3  SRR35899${i} &;done这个命令报错~~~。好像漏了后缀。(我从网上弄了一个别的)
        for i in *sra
        do
        echo $i
        fastq-dump --split-3 $i
        done
回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-10-14 22:24 , Processed in 0.042234 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.