宏转录组学习笔记（⼆）
上次竟然忘记把学习的教程地址放上了，这⾥⾸先放上：
继续前⾯的学习，前⾯已经把软件安装完成，数据库准备好，下⾯就是分析的过程了，基本上按照原⽂的命令进⾏的，由于教程中没有给出tara_f135_full_megahit.fasta这个⽂件，这⾥我们就把这⼏个样本的组装提到了前⾯，⾃⼰组装获得这个序列，然后再进⾏物种注释。
质控
把不合格的数据去除，这⾥使⽤了fastqc对每个样本进⾏检查，multiqc汇兑结果，然后Trimmomatic进⾏去除。
1.fastqc检查，multiqc汇兑结果
#激活⼯作环境
conda activate tera
#为了⽅便，设置变量，临时有⽤，退出失效，所以再运⾏⼀次
过生日送什么礼物
export PROJECT=~/work
#进⼊⽂件夹
cd${PROJECT}
#建⽴质控⽂件夹并进⼊，每个步骤建⽴⽂件夹是个好习惯
mkdir -p quality
cd quality
#软链接数据，省去路径问题，也是⼀个好⽅法
ln -s ../data/*.fq.gz ./
#再次check⼀切就绪
printf"I see $(ls -1 *.fq.gz |wc -l) files here.\n"
#开始质控，原⽂没有加⼊多线程，这⾥我使⽤4线程
fastqc *.fq.gz -t 4
#查看确认运⾏完成，⼀般没有问题
ls -d *fastqc.zip*
'''
TARA_135_DCM_5-20_rep1_1m_1_fastqc.zip TARA_135_SRF_5-20_rep2_1m_1_fastqc.zip
...
TARA_135_SRF_5-20_rep1_1m_2_fastqc.zip TARA_136_SRF_5-20_rep2_1m_2_fastqc.zip
'''
#multiqc汇兑结果
multiqc .
运⾏完成后使⽤浏览器查看下结果如TARA_135_SRF_5-20_rep1_1m_1_fastqc.html，就是我们熟悉的⽹页了。
2.根据质量修剪和过滤
#建⽴新的⽂件夹，软链⽂件，准备adapters
cd$PROJECT
mkdir -p trim
cd trim
ln -s ../data/*.fq.gz .
cat ~/Miniconda/envs/tera/share/trimmomatic-*/adapters/* > combined.fa
#⼀个⼤的for循环解决问题
for filename in *;
do
#使⽤ basename 变量去除 _ 产⽣ base变量，第⼀次见这种操作
base=$(basename $filename _)
echo$base
# 运⾏ Trimmomatic
trimmomatic PE ${base}_ \
${base}_ \
${base}_1. s1_se \
${base}_2. s2_se \
ILLUMINACLIP:combined.fa:2:40:15 \
LEADING:2 TRAILING:2 \
SLIDINGWINDOW:4:2 \
MINLEN:25
# 保存orphans序列
gzip -9c s1_se s2_se >> orphans.
rm -f s1_se s2_se
done
#看下质控效果，同样是fastqc-multiqc，⼀起完成，稍微看了眼，只去掉了很少序列数
fastqc *. -t 4
multiqc .
#取消⽬录下的所有⽂件可写权限
chmod a-w ${PROJECT}/trim/*.
更多质控trim策略看这， ，这是转录组的教程，但是原理应该通⽤。
Khmer错误修剪
前⾯的trimmomatic修剪后，依然有序列包含错误，Q30意味着1000个Q值是30的碱基中，约有⼀个是错的。⾼质量碱基也会⼩概率的错误，反过来，许多低质量碱基也是正确的。然后，我们修剪地很轻微，只去除了质量极低的碱基。这是有意为之，因为组装需要保留尽可能⼤的测序深度，组装软件可以通过测序深度判断正确与否。另⼀个修剪选择是基于k-mer丰度（kmer spectral error），⽂章显⽰结果优于基于质量修剪。
它的逻辑是这样的：如果你在⾼测序深度发现低丰度的K-mers，这些很可能就是错误。当然，菌株变异也会引起这样。在宏转录组中可能遇到极低丰度的菌的情况，所以我们不必须做这⼀步。我们实验室开发了⼀个⽅法，把序列按丰度⾼低分成两个组，仅修剪⾼丰度的序列，这样可以兼得低丰度序列，⼜能去除错误。
#建⽴⽂件夹并进⼊
cd${PROJECT}
mkdir -p khmer_trim
cd khmer_trim
#选择⼀个样本开始，链接⽂件
ln -s ${PROJECT}/trim/TARA_135_SRF_5-20_* ./
#还是for循环运⾏，这⼀步是最耗内存和时间的，运⾏需要8G以上内存，否则会使⽤swap，运⾏速度变得超慢
#实测在8G内存上，如果有图形界⾯的系统会卡死，开的⽆图形界⾯的去服务器使⽤内存7.8G
#参数⾥已经设置了内存8G，可能改⼩点也能运⾏，速度差别，没有测试
for filename in *_1.
do
#Use the program basename to remove _1. to generate the base
base=$(basename $filename _1.)
echo$base健康茶
interleave-reads.py ${base}_1. ${base}_2. | \
trim-low-abund.py - -V -Z 10 -C 3 -o - --gzip -M 8e9 | \
extract-paired-reads.py --gzip -p ${base}.khmer. -s ${base}.khmer.
done
'''
removing ./tmph8zrazrjkhmer/-.pass2
removing temp directory & contents (./tmph8zrazrjkhmer)
read 1995938 reads, 194677930 bp
wrote 1980263 reads, 192516464 bp
looked at 1910590 reads twice (1.96 passes)
removed 15675 reads and trimmed 36041 reads (2.59%)
高手如云trimmed or removed 1.11%% of bases (2161466 total)
1995938 reads were high coverage (100.00%);
skipped 0 reads/0 bases because of low coverage
fp rate estimated to be 0.000
output streamed to stdout
DONE; read 1980263 sequences, 982692 pairs and 14879 singletons
wrote to: TARA_135_SRF_5-20_rep2_1m.khmer. and TARA_135_SRF_5-20_rep2_1m.khmer.
'''
#查看k-mer丰度变化，去除了少部分序列（⼏万条）
unique-kmers.py TARA_135_SRF_5-20_rep1_1m_1. TARA_135_SRF_5-20_rep1_1m_2.
'''
Estimated number of unique 32-mers in TARA_135_SRF_5-20_rep1_1m_1.: 58871706
Estimated number of unique 32-mers in TARA_135_SRF_5-20_rep1_1m_2.: 57626073
Total estimated number of unique 32-mers: 100687408
'''
unique-kmers.py TARA_135_SRF_5-20_rep1_1m.khmer.
'''
Estimated number of unique 32-mers in TARA_135_SRF_5-20_rep1_1m.khmer.: 100296061
Total estimated number of unique 32-mers: 100296061
赵丽颖陷害袁姗姗'''
因为这是专门⽤来培训的⼩数据，如果真实数据去除的序列数应该⽐这个⼤得多。有更多相关知识。
宏转录组组装
基本过程就是你丢给软件⼀些reads，获得⼀些代表你的转录组的contigs，或者拉长没有太多明显多态性和长重复的reads。你使⽤修剪好的reads运⾏⼀个转录组组装程序，得到⼀些组装好的RNA。这些contigs代表环境中真核⽣物中的转录本（Poly-A mRNA）。
#建⽴⽂件夹，链接⽂件
cd$PROJECT
mkdir -p assembly
cd assembly
ln -fs ${PROJECT}/khmer_trim/*.khmer. .
ls
#组装，这⾥我采⽤12G内存，2线程，已经是我的笔记本的极限，实际运⾏下来只使⽤了⼏百M的样⼦
time megahit --12 TARA_135_SRF_5-20_rep1_1m.khmer.,TARA_135_SRF_5-20_rep2_1m.
khmer. --memory 12e9 --num-cpu-threads 4 --out-prefix TARA_135_SRF --out-dir ./TARA_135_SRF_khmer -f
详细输出如下：
'''
2020-03-15 13:29:51 - MEGAHIT v1.2.9
2020-03-15 13:29:51 - Using megahit_core with POPCNT and BMI2 support
2020-03-15 13:29:51 - Convert reads to binary library
2020-03-15 13:29:55 - b'INFO sequence/io/sequence_lib.cpp : 77 - Lib 0 (/home/ubuntu/work/assembly/TARA_135_SRF_5-20_rep1_1m.khmer.pe.fq.g z): interleaved, 1958670 reads, 101 max length'
2020-03-15 13:29:58 - b'INFO sequence/io/sequence_lib.cpp : 77 - Lib 1 (/home/ubuntu/work/assembly/TARA_135_SRF_5-20_rep2_1m.khmer.pe.fq.g z): interleaved, 1965384 reads, 101 max length'
2020-03-15 13:29:58 - b'INFO utils/utils.h : 152 - Real: 7.6951\tuser: 2.6234\tsys: 0.4533\tmaxrss: 102644'
2020-03-15 13:29:58 - k-max reset to: 119
2020-03-15 13:29:58 - Start assembly. Number of CPU threads 4
2020-03-15 13:29:58 - k list: 21,29,39,59,79,99,119
2020-03-15 13:29:58 - Memory used: 8000000000
2020-03-15 13:29:58 - Extract solid (k+1)-mers for k = 21
2020-03-15 13:30:55 - Build graph for k = 21
......
2020-03-15 13:42:29 - Build graph for k = 119
2020-03-15 13:42:31 - Assemble contigs from SdBG for k = 119
2020-03-15 13:42:41 - Merging to output final contigs
2020-03-15 13:42:41 - 16707 contigs, total 8656583 bp, min 202 bp, max 9085 bp, avg 518 bp, N50 520 bp
2020-03-15 13:42:41 - ALL DONE. Time elapsed: 770.087650 seconds
奶粉十大品牌排行榜#时间有点长，有⼗⼆分钟。可以休息下，或者继续阅读后⾯的内容，准备下。
real 12m50.137s
user 23m11.702s
sys 0m7.025s
'''
q q 邮箱登陆查看组装⽂件
输出⽂件会是TARA_135_SRF_khmer/TARA_igs.fa。
#复制并重命名组装⽂件
cp ./TARA_135_SRF_khmer/*contigs.fa tara135_SRF_megahit.fasta
head tara135_SRF_megahit.fasta
我们为什么组装，组装的优点是什么，对于组装⽂件，我们能对它做些什么呢？
组装集减少读取冗余，因此组装每个转录本应具有⼤约⼀个序列，⽽不是读取，因为读取每个脚本具有许多读取。
组装的连续体⽐读取时间长，因此更容易对它们进⾏基因搜索。我们明天再报道！
程序集的错误也⽐读取少，因此样本⽐较等可能更准确。但是，如果数据覆盖范围确实很低，并且⼤量信息也会丢失，则程序集也可能消除某些数据。
进⼀步参考
还有其他的组装软件，你可以⽤您的数据尝试，其中⼀些列在我们的。我们感兴趣的⼀个是它在蛋⽩质⽔平组装！
组装评估和定量
是⼀个可⽤于⼏种不同类型的组装评估的程序。第⼀个也是最简单的⽅法是计算有关组装基于长度的指标，例如总碱基数和contigs的
N50。Transrate还可⽤于⽐较两个组装或给出⼀个分数，该分数表⽰为组装提供积极⽀持的输⼊reads的⽐例。有关指标以及如何运⾏基于参考的Transrate的更多信息，请参阅和论⽂。
#检查$PROJECT变量
echo$PROJECT
cd$PROJECT
#建⽴⽂件夹并进⼊
mkdir -p evaluation
cd evaluation
#把组装⽂件复制过来，应该是⽂件⽐较⼩就不链接了
cp$PROJECT/assembly/tara135_SRF_megahit.fasta ./
#使⽤transrate产⽣组装统计数据，这⾥我没有加环境变量，直接绝对路径运⾏了
./transrate-1.0.3-osx/transrate -h #查看帮助，测试软件可⽤
#统计
time ./transrate-1.0.3-osx/transrate --assembly tara135_SRF_megahit.fasta
'''
[ INFO] 2020-03-15 14:15:45 : Loading assembly: /Users/zd200572/work/evaluation/evaluation/tara135_SRF_megahit.fasta
[ INFO] 2020-03-15 14:15:46 : Analysing assembly: /Users/zd200572/work/evaluation/evaluation/tara135_SRF_megahit.fasta
[ INFO] 2020-03-15 14:15:46 : Results will be saved in /Users/zd200572/work/evaluation/evaluation/transrate_results/tara135_SRF_megahit
[ INFO] 2020-03-15 14:15:46 : Calculating
[ INFO] 2020-03-15 14:15:48 : Contig metrics:
[ INFO] 2020-03-15 14:15:48 : -----------------------------------
[ INFO] 2020-03-15 14:15:48 : n seqs 16707
[ INFO] 2020-03-15 14:15:48 : smallest 202
[ INFO] 2020-03-15 14:15:48 : largest 9085
[ INFO] 2020-03-15 14:15:48 : n bases 8656583
[ INFO] 2020-03-15 14:15:48 : mean len 518.14
[ INFO] 2020-03-15 14:15:48 : n under 200 0
[ INFO] 2020-03-15 14:15:48 : n over 1k 909
[ INFO] 2020-03-15 14:15:48 : n over 10k 0
[ INFO] 2020-03-15 14:15:48 : n with orf 5531
[ INFO] 2020-03-15 14:15:48 : mean orf percent 84.25
[ INFO] 2020-03-15 14:15:48 : n90 331
[ INFO] 2020-03-15 14:15:48 : n70 408
[ INFO] 2020-03-15 14:15:48 : n50 520
[ INFO] 2020-03-15 14:15:48 : n30 711
[ INFO] 2020-03-15 14:15:48 : n10 1208
[ INFO] 2020-03-15 14:15:48 : gc 0.5
[ INFO] 2020-03-15 14:15:48 : bases n 0
[ INFO] 2020-03-15 14:15:48 : proportion n 0.0
[ INFO] 2020-03-15 14:15:48 : Contig metrics done in 2 seconds
[ INFO] 2020-03-15 14:15:48 : No reads provided, skipping read diagnostics
[ INFO] 2020-03-15 14:15:48 : No reference provided, skipping comparative diagnostics
[ INFO] 2020-03-15 14:15:48 : Writing contig metrics for each contig to /Users/zd200572/work/evaluation/evaluation/transrate_results/tara135_SRF_mega hit/contigs.csv
[ INFO] 2020-03-15 14:15:48 : Writing analysis results to assemblies.csv
real 0m4.218s
user 0m3.837s
sys 0m0.291s
'''
ln -s ../assembly/tara_f135_full_megahit.fasta
# full vs. subset 因为这⾥没有全部的，这步没做
./transrate-1.0.3-osx/transrate --reference=tara_f135_full_megahit.fasta --assembly=tara135_SRF_megahit.fasta --output=full_v_subset
# subset vs. full 因为这⾥没有全部的，这步没做
./transrate-1.0.3-osx/transrate --assembly=tara_f135_full_megahit.fasta --reference=tara135_SRF_megahit.fasta --output=subset_v_full
⽤salmon量化reads
这个组装对我们测序reads的代表性如何，我们将使⽤来量化表达。Salmon 是⼀种⽤于量化 RNAseq reads的新软件，它既运⾏快，⼜考虑了转录本的长度）。
您可以在此处阅读更多有关salmon类似的"pseudoalignment"：
博客⽂章简介： robpatro/blog/?p=248
2016年博客⽂章评估和⽐较
salmon
2018年论⽂：[](/content/early/2018/10/16/444620)
#⾸先链接⽂件
ln -s ${PROJECT}/trim/TARA_135_SRF_5-20_*. ./
#⽤salmon进⾏定量
版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。
发表评论