八 09

RNAseq数据完整生物信息分析流程第一讲之文献数据下载

Posted on 2016年8月9日 by ulwvfje

我这里拿的是bioconductor里面最常用的airway数据，因为差异表达分析在bioconductor里面是重点，它们这些包在介绍自己的算法以及做示范的时候都用的这个数据。可以在GEO数据库里面看到信息描述：http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52778 可以看到是Illumina HiSeq 2000 (Homo sapiens) ，75bp paired-end 这个信息很重要，决定了下载sra数据之后如何解压以及如何比对。也可以看到作者把所有的测序原始数据都上传到了SRA中心：http://www.ncbi.nlm.nih.gov/sra?term=SRP033351 ，这里可以在linux服务器上面写一个简单的脚本批量下载所有的测序数据，然后根据GEO里面描述的metadata把原始数据改名。

for ((i=508;i<=523;i++)) ;do wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP033/SRP033351/SRR1039$i/SRR1039$i.sra;done
ls *sra |while read id; do ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 $id;done

需要自己看SRA里面的数据记录，上面的脚本不难写出，然后因为是Illumina的双端测序，所以我们用fastq-dump --split-3命令来把sra格式数据转换为fastq，但是因为这里有16个测序数据，所以最好是同步改名，我这里用脚本批量生成改名脚本如下：

为了节省空间，我用了--gzip压缩，该文件名，用-A参数。

nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N61311_untreated SRR1039508.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N61311_Dex SRR1039509.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N61311_Alb SRR1039510.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N61311_Alb_Dex SRR1039511.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N052611_untreated SRR1039512.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N052611_Dex SRR1039513.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N052611_Alb SRR1039514.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N052611_Alb_Dex SRR1039515.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N080611_untreated SRR1039516.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N080611_Dex SRR1039517.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N080611_Alb SRR1039518.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N080611_Alb_Dex SRR1039519.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N061011_untreated SRR1039520.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N061011_Dex SRR1039521.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N061011_Alb SRR1039522.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N061011_Alb_Dex SRR1039523.sra &

可以看到这里的16个样本来源于同样的4个人，是HASM细胞系，处理详情如下：

测序基础：

HASM细胞系-human airway smooth muscle，

The Illumina TruSeq assay was used to prepare 75bp paired-end libraries for HASM cells from four white male donors under four treatment conditions:

1) no treatment;

2) treatment with a β2-agonist (i.e. Albuterol, 1μM for 18h);

3) treatment with a glucocorticosteroid (i.e. Dexamethasone (Dex), 1μM for 18h);

4) simultaneous treatment with a β2-agonist and glucocorticoid

and the libraries were sequenced with an Illumina Hi-Seq 2000 instrument.

我们这里只是先根据fastq数据比对到参考基因组，然后计算每个样本的表达量即可，后续的分组计算差异表达，就需要个性化了。

下载的sra大小如下：

-rw-rw-r-- 1 jmzeng jmzeng 1.6G Aug 9 04:21 SRR1039508.sra
-rw-rw-r-- 1 jmzeng jmzeng 1.5G Aug 9 05:20 SRR1039509.sra
-rw-rw-r-- 1 jmzeng jmzeng 1.6G Aug 9 06:14 SRR1039510.sra
-rw-rw-r-- 1 jmzeng jmzeng 1.5G Aug 9 07:05 SRR1039511.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.1G Aug 9 08:07 SRR1039512.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.3G Aug 9 09:17 SRR1039513.sra
-rw-rw-r-- 1 jmzeng jmzeng 3.1G Aug 9 10:56 SRR1039514.sra
-rw-rw-r-- 1 jmzeng jmzeng 1.9G Aug 9 11:56 SRR1039515.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.1G Aug 9 13:02 SRR1039516.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.6G Aug 9 14:16 SRR1039517.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.3G Aug 9 15:17 SRR1039518.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.0G Aug 9 16:05 SRR1039519.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.1G Aug 9 16:56 SRR1039520.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.4G Aug 9 17:57 SRR1039521.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.0G Aug 9 18:46 SRR1039522.sra
-rw-rw-r-- 1 jmzeng jmzeng 1.4G Aug 9 19:28 SRR1039523.sra

解压后成双端测序的fastq数据如下：

-rw-rw-r-- 1 jmzeng jmzeng 2.5G Aug 9 20:12 N052611_Alb_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 2.5G Aug 9 20:12 N052611_Alb_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 20:44 N052611_Alb_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 20:44 N052611_Alb_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 289M Aug 9 20:44 N052611_Alb_Dex.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 951M Aug 9 20:59 N052611_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 954M Aug 9 20:59 N052611_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.7G Aug 9 20:53 N052611_untreated_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.7G Aug 9 20:53 N052611_untreated_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.5G Aug 9 20:45 N061011_Alb_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.5G Aug 9 20:45 N061011_Alb_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.9G Aug 9 20:59 N061011_Alb_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.9G Aug 9 20:59 N061011_Alb_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 16M Aug 9 20:45 N061011_Alb.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.4G Aug 9 20:48 N061011_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.4G Aug 9 20:48 N061011_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.2G Aug 9 20:00 N061011_untreated_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.2G Aug 9 20:00 N061011_untreated_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 759M Aug 9 20:00 N061011_untreated.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.9G Aug 9 20:03 N080611_Alb_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.9G Aug 9 20:03 N080611_Alb_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 19:59 N080611_Alb_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 19:59 N080611_Alb_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 535M Aug 9 19:59 N080611_Alb_Dex.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 2.1G Aug 9 20:06 N080611_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 2.1G Aug 9 20:06 N080611_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.6G Aug 9 20:01 N080611_untreated_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.6G Aug 9 20:01 N080611_untreated_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:09 N61311_Alb_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:09 N61311_Alb_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:08 N61311_Alb_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:08 N61311_Alb_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.2G Aug 9 08:07 N61311_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.2G Aug 9 08:07 N61311_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:09 N61311_untreated_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:09 N61311_untreated_2.fastq.gz

接下来所有的分析就基于此数据啦

七 05

自学CHIP-seq分析第三讲~公共测序数据下载

Posted on 2016年7月5日 by ulwvfje

这一步跟自学其它高通量测序数据处理一样，就是仔细研读paper，在里面找到作者把原始测序数据放在了哪个公共数据库里面，一般是NCBI的GEO，SRA，本文也不例外，然后解析样本数，找到下载链接规律

## step1 : download raw data

cd ~

mkdir CHIPseq_test && cd CHIPseq_test

mkdir rawData && cd rawData

## batch download the raw data by shell script :

for ((i=593;i<601;i++)) ;do wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP033/SRP033492/SRR1042$i/SRR1042$i.sra;done

很容易就下载了8个测序文件，每个样本的数据大小，测序量如下

621M Jun 27 14:03 SRR1042593.sra (16.9M reads)

2.2G Jun 27 15:58 SRR1042594.sra (60.6M reads)

541M Jun 27 16:26 SRR1042595.sra (14.6M reads)

2.4G Jun 27 18:24 SRR1042596.sra (65.9M reads)

814M Jun 27 18:59 SRR1042597.sra (22.2M reads)

2.1G Jun 27 20:30 SRR1042598.sra (58.1M reads)

883M Jun 27 21:08 SRR1042599.sra (24.0M reads)

2.8G Jun 28 11:53 SRR1042600.sra (76.4M reads)

虽然下载的SRA格式数据也是一个很流行的标准，但它只是数据压缩的标准，几乎没有软件能直接跟SRA的格式的测序数据来进行分析，我们需要转成fastq格式，代码如下：

## step2 : change sra data to fastq files.

## cell line: MCF7 // Illumina HiSeq 2000 // 50bp // Single ends // phred+33

## http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52964

## ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP033/SRP033492

ls *sra |while read id; do ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump $id;done

rm *sra

解压的详情如下，可以看到SRA格式有6~9倍的压缩了，比zip格式压缩的2~3倍高多了

## 621M --> 3.9G

## 2.2G --> 14G

## 541M --> 3.3G

## 2.4G --> 15G

三 26

阅读捕获测序文章并下载数据

Posted on 2015年3月26日 by ulwvfje

一．阅读文献找到SRP

该文献讲了单分子测序在医疗领域的一个应用，我感觉挺重要的，就分析了一下，然后下载了数据，准备处理一下。

Single-step capture and sequencing of natural DNA for detection of BRCA1 mutations

在NCBI查到该数据地址，并且用脚本下载即可

http://www.ncbi.nlm.nih.gov/sra/?term=SRP007097

下载之后的数据如下，共19个测序文件，都是200K左右大小，那两个一百多M的可能是下载错了

for i in {32..52}

wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR258/SRR2588$i/SRR2588$i.sra

Done

下载的19个数据，都是只有1万多条序列。

因为这些判断都是对BRCA1这个基因进行目标性测序，所以接下来需要对它们进行特殊的处理。

三 17

草莓基因组文章解读-并下载原始测序数据

Posted on 2015年3月17日 by ulwvfje

找橡胶测序数据无果

所以我只好找了他们所参考的草莓（strawberry, Fragaria vesca (2n = 2x = 14)，a small genome (240 Mb),）的文章，是发表是nature genetics上面的

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3326587/

可以看到它的SRA索取号。

草莓组装结果：Over 3,200 scaffolds were assembled with an N50 of 1.3 Mb .

Over 95% (209.8 Mb) of the total sequence is represented in 272 scaffolds.

草莓基因息：Gene prediction modeling identified 34,809 genes, with most being supported by transcriptome mapping.

草莓染色体信息：Paradoxically, the small basic (x = 7) genome size of the strawberry genus, ~240 Mb,

offers substantial advantages for genomic research.

草莓来源：diploid strawberry F. vesca ssp. vesca accession Hawaii 4

(National Clonal Germplasm Repository accession # PI551572).

然后我去NCBI上面下载这三个数据

SRA020125 共有四个数据：

http://www.ncbi.nlm.nih.gov/sra/SRX030575[accn]	Total: 4 runs, 4.7M spots, 2.6G bases, 5.5Gb
http://www.ncbi.nlm.nih.gov/sra/SRX030576[accn] （3 KB PE）	Total: 2 runs, 2.2M spots, 908.5M bases, 2.1Gb
http://www.ncbi.nlm.nih.gov/sra/SRX030577[accn] （20KB片段）	Total: 2 runs, 1.9M spots, 800M bases, 1.8Gb
http://www.ncbi.nlm.nih.gov/sra/SRX030578[accn]	Total: 3 runs, 4M spots, 2.2G bases, 4.6Gb

挂在后台自动下载

好了，有了这些数据我们就要进行基因组的一系列分析啦！！！

不过我们可以先看看他们这个研究小组的成果

首先他们建造了一个关于草莓的基因组信息网站

https://strawberry.plantandfood.co.nz/

跟我之前在水科院做鲫鱼鲤鱼的差不多

直接在里面就可以下载他们做好的所有数据，也可以可视化。

它的染色体如下，非常简单，就七条染色体

http://www.rosaceae.org/species/fragaria/fragaria_vesca/genome_v1.1

我找到了它组装好的草莓基因组地址，用批处理全部下载了

三 17

研读橡胶的基因组文章-结果没有原始测序数据

Posted on 2015年3月17日 by ulwvfje

研读橡胶的基因组文章

我本科的前两年在海南儋州读书，那时候旁边就是橡胶所，很多同学也在那边做毕业论文什么的，我一直以为那里是全世界的橡胶中心，所有的先进技术都在那里产生，结果，前些天跟一个橡胶所的老师聊天才发现，居然橡胶(Hevea brasiliensis)的基因组已经发表了，可是，跟橡胶所没有半毛钱关系，更搞笑的事情是，堂堂一个基因组文章居然发表在BMC这样的杂志，真不知道是基因组的年代已经过去了还是他们做的实在是太差了，反正我看不过去了，所以研读他们的文章，并且下载数据测试一下。

文章地址如下:http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3575267/

可以看到它过于数据的描述都在补充材料1里面，所以我下载了补充材料。

可以看到所有的测序数据的描述，45个G的i llumina的200bp的双端测序，27个G的illumina的200bp的双端测序，约10G左右的长片段（8kb，20kb）罗氏454数据，最后还有一点点solid数据，它这样的测序策略好像是模仿的2011年发布的草莓基因组数据。

但是补充材料里面没有列出下载地址，我有点困惑！

按照道理我研读文献的步骤应该没有错，有可能是因为这个文章发表的杂志水平太低，所以不要求他们把测序原始数据上传到NCBI的SRA里面。或者是他们本身觉得文章发的不够档次，不想公布数据，所以先留着自己做精细分析，等发了大文章再公布原始数据。

然后我在NCBI的SRA里面查找了关于橡胶的原始数据，果真没有

仅有的10个数据，都是别的小组做的RNA-seq的内容。

De novo transcriptome analysis of abiotic stress responsive transcripts of Hevea brasiliensis.