无参转录组起步

  今天是生信星球陪你的第266天

---

   大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~

   就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~

   这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！

刘小泽写于19.2.1

😂发觉好久没写了，果然临近过年事情比较多，最近重新开始无参转录组。

上次写了个开篇，这次先从数据下载做到得到组装的转录本

数据来源

Defining the transcriptomic landscape of Candida glabrata by RNA-Seq.

http://www.ncbi.nlm.nih.gov/pubmed/?term=25586221

这里只下载GSNO和wt的数据SRR1582646-SRR1582651便于分析

需要注意的是，如果用自己的数据没有问题，要是使用公共数据的sra再转换，就需要注意一个问题：不要使用平时默认的参数--gzip --split-3
这样会产生这种header：@SRR1582646.224153461 ，而这种是不被软件识别的

我们要使用参数：--defline-seq '@$sn[_$rn]/$ri' --split-files ，这样得到正确的header结果：@HWI-ST1363:114:D1H5CACXX:3:1101:1196:2196/1

数据下载转换

SRR_Acc_List.txt中放入SRR的ID号

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 [e-mail]:/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 


预处理（质控+过滤）

因为是测试数据，因此我们可以只选择其中10k行，但是为了选取的随机性，可以使用seqtk软件

source activate rna-seq
# 选择一部分10000条，用seqtk sample随机选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

第一次质控

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

再一次质控

cd $wkd/clean
source activate rna-seq
fastqc -o ./ -t 10 $wkd/clean/*.fq
multiqc ./
source deactivate


Trinity拼接

如果分析的数据多，就可以做个表格；数据少的话，直接指定--left，--right

这里在正式分析前，先做一个样本表格samples.txt，其中包含了处理、重复、左端数据、右端数据

# 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


跑程序很简单，需要指定样本文件、CPU、内存、min_contig_length(minimum assembled contig length)【指定拼接最小长度，默认200bp，低于这个长度的contig就被丢弃】

Trinity --seqType fq  --samples_file $wkd/assembly/samples.txt \
      --CPU 10 --max_memory 10G --min_contig_length 150

运行Trinity中间会用到samtools，但是有时候安装它会失败，会遇到动态库的问题（就是某些lib缺失）。最简单的办法就是降级，如最新版1.9，我们安装1.8 ，例如conda install samtools=1.8

剩下的运行过程中 ~1 hour and ~1G RAM per ~1 million PE  reads

首先运行第一个程序：seqtk-trinity ，目的是将fq转为fa，并生成both.fa

第二个程序：Jellyfish ，中间过程jellyfish count + histo + dump，starting with Jellyfish to generate the k-mer catalog

第三个程序：Inchworm， assemble 'draft' contigs

第四个程序：Chrysalis，cluster the contigs and build de Bruijn graphs

第五个程序：Butterfly，tracing paths through the graphs and reconstructing the final isoform sequences

最后的结果就是得到了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信息中包含了：转录本名称 、长度、 de Bruijn图重构路径，其中g表示gene，i表示isoform转录本，多个转录本可以重构成一个基因，在下游差异分析中，可以基于基因或转录本两种层次进行分析

---