劉小澤寫于19.2.1
??發(fā)覺好久沒寫了，果然臨近過年事情比較多鉴逞，最近重新開始無參轉(zhuǎn)錄組。這次先從數(shù)據(jù)下載做到得到組裝的轉(zhuǎn)錄本
數(shù)據(jù)來源
Defining the transcriptomic landscape of Candida glabrata by RNA-Seq.
http://www.ncbi.nlm.nih.gov/pubmed/?term=25586221
這里只下載GSNO和wt的數(shù)據(jù)SRR1582646-SRR1582651
需要注意的是联喘，如果用自己的數(shù)據(jù)沒有問題华蜒，要是使用公共數(shù)據(jù)的sra再轉(zhuǎn)換辙纬，就需要注意一個(gè)問題：不要使用平時(shí)默認(rèn)的參數(shù)
--gzip --split-3
這樣會(huì)產(chǎn)生這種header：@SRR1582646.224153461豁遭，而這種是不被軟件識(shí)別的我們要使用參數(shù)：
--defline-seq '@$sn[_$rn]/$ri' --split-files，這樣得到正確的header結(jié)果：@HWI-ST1363:114:D1H5CACXX:3:1101:1196:2196/1
數(shù)據(jù)下載轉(zhuǎn)換
SRR_Acc_List.txt中放入SRR的ID號(hào)
SRR1582646
SRR1582647
SRR1582648
SRR1582649
SRR1582650
SRR1582651
wkd=~/rna-denovo-test
source activate rna-seq
mkdir $wkd/raw && cd $wkd/raw
cat $wkd/raw/SRR_Acc_List.txt |while read i
do
ascp -v -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh -k 1 -T \
-l200m anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/SRR/SRR158/$i/${i}.sra ./
time fastq-dump --defline-seq '@$sn[_$rn]/$ri' --split-files -A $i ${i}.sra && echo "** ${i}.sra to fastq done **"
done
source deactivate
預(yù)處理（質(zhì)控+過濾）
因?yàn)槭菧y(cè)試數(shù)據(jù)贺拣，因此我們可以只選擇其中10k行蓖谢，但是為了選取的隨機(jī)性，可以使用seqtk軟件
source activate rna-seq
# 選擇一部分10000條譬涡，用seqtk sample隨機(jī)選10000條
ls $wkd/raw/*_1.fastq.gz >1
ls $wkd/raw/*_2.fastq.gz >2
paste 1 2 > sample-conf
cat conf1 | while read i;do
fqs=($i)
fq1=${fqs[1]}
fq2=${fqs[2]}
surname=${fq1%%_*}
#echo $fq1 $surname
seqtk sample $fq1 10000 >test.${surname}_1.fastq
seqtk sample $fq2 10000 >test.${surname}_2.fastq
done
第一次質(zhì)控
mkdir $wkd/qc && cd $wkd/qc
fastqc -o ./ -t 10 $wkd/raw/*.fastq
multiqc ./
source deactivate
過濾
source activate rna-seq
mkdir $wkd/clean && cd $wkd/clean
ls $wkd/raw/*_1.fastq >1
ls $wkd/raw/*_2.fastq >2
paste 1 2 > filter_conf
cat filter_conf | while read i;do
fqs=($i)
fq1=${fqs[0]}
fq2=${fqs[1]}
trim_galore -q 25 --phred33 --length 50 -e 0.1 --stringency 3 \
--paired -o $wkd/clean $fq1 $fq2
done
source deactivate
再一次質(zhì)控
cd $wkd/clean
source activate rna-seq
fastqc -o ./ -t 10 $wkd/clean/*.fq
multiqc ./
source deactivate
Trinity拼接
如果分析的數(shù)據(jù)多闪幽，就可以做個(gè)表格；數(shù)據(jù)少的話涡匀，直接指定--left盯腌，--right
這里在正式分析前，先做一個(gè)樣本表格samples.txt陨瘩，其中包含了處理腕够、重復(fù)级乍、左端數(shù)據(jù)、右端數(shù)據(jù)
# cond_A cond_A_rep1 A_rep1_left.fq A_rep1_right.fq
# cond_A cond_A_rep2 A_rep2_left.fq A_rep2_right.fq
# cond_B cond_B_rep1 B_rep1_left.fq B_rep1_right.fq
# cond_B cond_B_rep2 B_rep2_left.fq B_rep2_right.fq
mkdir $wkd/assembly && cd $wkd/assembly
echo -e "GSNO\nGSNO\nGSNO\nwt\nwt\nwt" >1
echo -e "GSNO_SRR1582646\nGSNO_SRR1582647\nGSNO_SRR1582648\nwt_SRR1582649\nwt_SRR1582650\nwt_SRR1582651" >2
ls $wkd/clean/*1_val_1.fq >3
ls $wkd/clean/*2_val_2.fq >4
paste 1 2 3 4 >samples.txt
跑程序很簡(jiǎn)單帚湘，需要指定樣本文件玫荣、CPU、內(nèi)存大诸、min_contig_length(minimum assembled contig length)【指定拼接最小長(zhǎng)度捅厂，默認(rèn)200bp，低于這個(gè)長(zhǎng)度的contig就被丟棄】
Trinity --seqType fq --samples_file $wkd/assembly/samples.txt \
--CPU 10 --max_memory 10G --min_contig_length 150
運(yùn)行Trinity中間會(huì)用到samtools资柔，但是有時(shí)候安裝它會(huì)失敗焙贷，會(huì)遇到動(dòng)態(tài)庫的問題（就是某些lib缺失）。最簡(jiǎn)單的辦法就是降級(jí)贿堰，如最新版1.9盈厘，我們安裝1.8 ，例如
conda install samtools=1.8
剩下的運(yùn)行過程中 ~1 hour and ~1G RAM per ~1 million PE reads
首先運(yùn)行第一個(gè)程序：seqtk-trinity 官边，目的是將fq轉(zhuǎn)為fa沸手，并生成both.fa
第二個(gè)程序：Jellyfish ，中間過程jellyfish count + histo + dump注簿，starting with Jellyfish to generate the k-mer catalog
第三個(gè)程序：Inchworm契吉， assemble 'draft' contigs
第四個(gè)程序：Chrysalis，cluster the contigs and build de Bruijn graphs
第五個(gè)程序：Butterfly诡渴，tracing paths through the graphs and reconstructing the final isoform sequences
最后的結(jié)果就是得到了Trinity.fasta捐晶，如
>TRINITY_DN506_c0_g1_i1 len=171 path=[149:0-170] [-1, 149, -2]
TGAGTATGGTTTTGCCGGTTTGGCTGTTGGTGCAGCTTTGAAGGGCCTAAAGCCAATTGT
TGAATTCATGTCATTCAACTTCTCCATGCAAGCCATTGACCATGTCGTTAACTCGGCAGC
AAAGACACATTATATGTCTGGTGGTACCCAAAAATGTCAAATCGTGTTCAG
>TRINITY_DN512_c0_g1_i1 len=168 path=[291:0-167] [-1, 291, -2]
ATATCAGCATTAGACAAAAGATTGTAAAGGATGGCATTAGGTGGTCGAAGTTTCAGGTCT
AAGAAACAGCAACTAGCATATGACAGGAGTTTTGCAGGCCGGTATCAGAAATTGCTGAGT
AAGAACCCATTCATATTCTTTGGACTCCCGTTTTGTGGAATGGTGGTG
其中header信息中包含了：轉(zhuǎn)錄本名稱 、長(zhǎng)度妄辩、 de Bruijn圖重構(gòu)路徑惑灵，其中g(shù)表示gene，i表示isoform轉(zhuǎn)錄本眼耀，多個(gè)轉(zhuǎn)錄本可以重構(gòu)成一個(gè)基因英支，在下游差異分析中，可以基于基因或轉(zhuǎn)錄本兩種層次進(jìn)行分析
歡迎關(guān)注我們的公眾號(hào)～_～
我們是兩個(gè)農(nóng)轉(zhuǎn)生信的小碩哮伟，打造生信星球干花，想讓它成為一個(gè)不拽術(shù)語、通俗易懂的生信知識(shí)平臺(tái)楞黄。需要幫助或提出意見請(qǐng)后臺(tái)留言或發(fā)送郵件到Bioplanet520@outlook.com
