rsem对转录本进行定量

最近在研究转录本,现在在下载数据,想起来自己有一个博客,就暂且来这里更新一下内容。

要想对转录本进行定量,首先需要下载它的转录组数据,将别人上传的SRR文件的名字整理在wheat.txt中,引用

prefetch --option-file wheat.txt

下载后通过sratoolkits将sra数据转化成fq格式

fastq-dump --split-3 /home/SRR3134001.sra --gzip -O /home/data

--split-3将sra数据的双端拆分成两个文件,原来单端并不会保存成两个文件。gzip将其压缩 -O为输出文件夹

才发现下载好慢啊,一会打算试一下fasterq-dump,据说非常快,应该没毛病。2

fasterq-dump -e 40 --split-3 /home/SRR3134001.sra --gzip -O /home/data

下载fq文件后,开始正式分析。
1、首先对转录组数据进行质控,这里运用fastp写了一个循环

for i in $(seq 1 2)
do
fastp -w 16 \
-i ../wheat${i}_1.fq.gz \
-I ../wheat${i}_2.fq.gz \
-o wheat${i}_clean_1.fq.gz \
-O wheat${i}_clean_2.fq.gz \
--html wheat${i}.html --json wheat${i}.json
done

2、在ensembl plant上下载小麦genome和gtf文件,小麦基因组也太大了,,吓到我了,,还好服务器不会说话,不然多少得骂我几句了。
在开始定量前,首先需要构建索引。

rsem-prepare-reference --gtf genome.gtf genome.fa reference_name -p 8

--gtf genome.gtf:输入基因组GTF注释文件。
genome.fa:基因组文件。
reference_name:索引名称。
-p:线程数。

构建索引后开始定量,我在rsem中直接调用star,在这里再写一个循环

for i in *_1.fq.gz
do
rsem-calculate-expression --paired-end -p 40 --star --star-gzipped-read-file ../02cleandata/${i} ../02cleandata/${i%_*}_2.fq.gz ./${i%_*}
done

--paired-end:表示输入的数据为双端测序数据。
这样就得到结果了。
genes.results和isoforms.results分别是基于基因和转录本水平的定量结果。
isoforms.results中包含了转录本ID,基因ID,转录本长度,有效长度,expected_count,TPM,FPKM和IsoPct(该转录本表达量占基因总表达量的百分比)。

上一篇:构建Zookeeper镜像


下一篇:偷懒方式安装cmake-3.14.5