rsem對轉錄本進行定量

最近在研究轉錄本，現在在下載數據，想起來自己有一個博客，就暫且來這里更新一下內容。

要想對轉錄本進行定量，首先需要下載它的轉錄組數據，將別人上傳的SRR文件的名字整理在wheat.txt中，引用

prefetch --option-file wheat.txt

下載后通過sratoolkits將sra數據轉化成fq格式

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

--split-3將sra數據的雙端拆分成兩個文件,原來單端並不會保存成兩個文件。gzip將其壓縮 -O為輸出文件夾

剛才用了一些fasterq-dump巨快，，

fasterq-dump --split-3 SRR1542405.sra -O /home/wheat

下載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（該轉錄本表達量占基因總表達量的百分比）。

更新：有一些大的基因組，star會提醒內存不夠，所以我又用bowtie2試了一下
構建索引：

rsem-prepare-reference -gtf ~/home/wheat/04test/wheat.gtf --bowtie2 ~/home/wheat/04test/wheat.fa ~/home/wheat/04test/wheat

運行rsem-calculate-expression(官方文檔寫的很全，我下載的是雙端測序數據)

rsem-calculate-expression --bowtie-path /home/miniconda3/envs/rna-seq/bin \
--fragment-length-mean 150.0 \
--fragment-length-sd 35.0 \
-p 20 \
--output-genome-bam \
--paired-end \
--calc-ci \
--ci-memory 1024 \
/home/wheat/04test/wheat1_1.fq \
/home/wheat/04test/wheat1_2.fq \
/home/wheat/04test/wheat \
wheat