组装策略
适用于设计多样本多物种的组装。例如100个样本，10个物种。这里如果想直接完成10个de nove组装，需要将所有样本数据放到一起后，通过样本信息表声明每个样本的物种名，再分别组装。因此涉及通过awk操作样本信息表。另外，多物种的分析中涉及到直系同源基因。此时通过unigene进行连接是最佳的。但在transdecoder预测ORF时，有时同一个转录本预测得到两个不同的ORF。所以最终是以预测得到的ORF为单元进行定量和分析。响应的CDS和PEP用以计算直系同源基因以及后续进化分析等。因此，这个过程中，从Trinity.fasta得到unigene，再以unigene为参考基因组，预测ORF得到gff3和gtf文件，进行有参组装。这种方式用来规避用Trinity.fasta为参考时，涉及到同一个基因得到不同注释，以及后续提取最长转录本时遇到的麻烦等。
awk操作样本信息表
样本信息表：
$ cat sample_list.txt
2023NJFGX01-SZ GX
2023NJZGX02-SZ HF
2023NJCYX05-SZ CY
2023NJZGX01-SZ HF
2023NJCYX303-SZ SP
2023NJCYX04-SZ CY
2023NJFGX04-SZ GX
2023NJBSX02-SZ HF
2023NJCYX205-SZ CY
2023NJZGX03-SZ CY
2023NJCYX302-SZ SP
按照第二列进行分配
aaa="SP"
cat sample_list.txt |
awk -v bbb=$aaa '{if($2 ==bbb) {print $1,$2}}'
2023NJCYX303-SZ SP
2023NJCYX302-SZ SP
事实上，我们想要将第一列的样本按文件夹分配
mkdir ${aaa}_dir;
cat sample_list.txt |
awk -v bbb=$aaa '{if($2 ==bbb) {print $1}}' |
while read line;
do
mv ${line}* ./${aaa}_dir;
done
将上述的aaa按照第二列的名称进行循环
cat sample_list.txt | awk '{print $2}' | sort | uniq |
while read line;
do
mkdir ${line}_dir;
cat sample_list.txt |
awk -v spname=$line '{if($2 ==spname) {print $1}}' |
while read line2;
do
mv ${line2}* ./${line}_dir;
done;
done
这样所有的样本被归到响应的文件夹里
Trinity组装（单个）
基本代码如下：
left_all=`echo *_1.clean.fq.gz`
right_all=`echo *_2.clean.fq.gz`
Trinity \
--seqType fq \
--max_memory 30G \
--left $left_all \
--right $right_all \
--CPU 12
会报错，反正经常这样报错：
Thread 2 terminated abnormally: Error, cmd: gunzip -c /mnt/d/Link_CQH/px22/CY_dir/2023NJCYX04-SZ_2.clean.fq.gz | fastool --illumina-trinity --to-fasta >> right.fa 2> /mnt/d/Link_CQH/px22/CY_dir/2023NJCYX04-SZ_2.clean.fq.gz.readcount died with ret 256 at /home/cqh/miniconda3/bin/Trinity line 2183.
类似的报错一般是先手动执行报错的这一步，获得left.fa和right.fa，后面再运行就对了。有时候会顺手直接cat到一起。
for fn in *1.clean.fq.gz
do
gunzip -c ./${fn} | fastool --illumina-trinity --to-fasta >> left.fa 2> ./${fn}.readcount
done
for fn in *2.clean.fq.gz
do
gunzip -c ./${fn} | fastool --illumina-trinity --to-fasta >> right.fa 2> ./${fn}.readcount
done
cat left.fa right.fa > both.fa
left_all=`echo *_1.clean.fq.gz`
right_all=`echo *_2.clean.fq.gz`
Trinity \
--seqType fq \
--max_memory 30G \
--left $left_all \
--right $right_all \
--CPU 12
获得unigene集
/home/cqh/miniconda3/opt/trinity-2.1.1/util/misc/get_longest_isoform_seq_per_trinity_gene.pl ./Trinity.fasta > unigene.fasta
修改名称
由于前面根据物种信息表，将物种名信息放在文件夹里。现提取文件夹的物种名信息，加入unigene的文件名文件内的序列中。
注意，“物种名|Trinity名”这种格式虽然很不错，但“|”在建立索引时会报错，所以舍弃了。
dirname=`pwd | awk 'BEGIN{FS="/"}{print $NF}' | awk 'BEGIN{FS="_"}{print $1}'`
cat unigene.fasta | awk '{print $1}' | sed "s/>/>${dirname}_/g" > ${dirname}_unigene.fasta
预测unigene集的开放阅读框，并获得gff3注释文件，通过gffread获得gtf文件。
TransDecoder.LongOrfs -t ${dirname}_unigene.fasta
cp ./${dirname}_unigene.fasta.transdecoder_dir/longest_orfs.gff3 ./${dirname}_unigene.gff3
gffread ${dirname}_unigene.gff3 -T -o ${dirname}_unigene.gtf
将unigene集以及gtf文件放到ref文件夹，以备有参组装。
mkdir ref
seqkit seq ${dirname}_unigene.fasta -w 1000 > ./ref/genome.fasta
cp ./${dirname}_unigene.gtf ./ref/genome.gtf
集合脚本
run_DIY_trinity.sh
我们首先有一个自动组装trinity的脚本，这个脚本仅完成到ORF注释，因为不确定每个流程中需要哪些步骤（未必每次都需要定量或者功能注释）
for fn in *1.clean.fq.gz
do
gunzip -c ./${fn} | fastool --illumina-trinity --to-fasta >> left.fa 2> ./${fn}.readcount
done
for fn in *2.clean.fq.gz
do
gunzip -c ./${fn} | fastool --illumina-trinity --to-fasta >> right.fa 2> ./${fn}.readcount
done
cat left.fa right.fa > both.fa
left_all=`echo *_1.clean.fq.gz`
right_all=`echo *_2.clean.fq.gz`
Trinity --seqType fq \
--max_memory 30G \
--left $left_all \
--right $right_all \
--CPU 12
# 并非重复，第一遍经常运行失败
Trinity --seqType fq \
--max_memory 30G \
--left $left_all \
--right $right_all \
--CPU 12
#unigene
/home/cqh/miniconda3/opt/trinity-2.1.1/util/misc/get_longest_isoform_seq_per_trinity_gene.pl ./trinity_out_dir/Trinity.fasta > ./trinity_out_dir/unigene.fasta
dirname=`pwd | awk 'BEGIN{FS="/"}{print $NF}' | awk 'BEGIN{FS="_"}{print $1}'`
cat ./trinity_out_dir/unigene.fasta | awk '{print $1}' | sed "s/>/>${dirname}_/g" > ${dirname}_unigene.fasta
#ORF注释
TransDecoder.LongOrfs -t ${dirname}_unigene.fasta
cp ./${dirname}_unigene.fasta.transdecoder_dir/longest_orfs.gff3 ./${dirname}_unigene.gff3
cp ./${dirname}_unigene.fasta.transdecoder_dir/longest_orfs.pep ./${dirname}_unigene.pep.fa
cp ./${dirname}_unigene.fasta.transdecoder_dir/longest_orfs.cds ./${dirname}_unigene.cds.fa
gffread ${dirname}_unigene.gff3 -T -o ${dirname}_unigene.gtf
run_DIY_subread.sh
有参组装脚本，实际上就是定量
#
mkdir ref
seqkit seq ${dirname}_unigene.fasta -w 1000 > ./ref/genome.fasta
cp ./${dirname}_unigene.gtf ./ref/genome.gtf
# index
path=`pwd`
subread-buildindex -o $path/ref/genome_index $path/ref/genome.fasta
# subjunc
index="$path/ref/genome_index"
for fn in *_1.clean.fq.gz
do
sample=${fn%_*}
subjunc -i $index \
-r ${sample}_1.clean.fq.gz \
-R ${sample}_2.clean.fq.gz \
-T 12 \
-o ${sample}.bam \
1>${sample}.subjunc.log 2>&1
done
#
mkdir clean_data
mv *.fq.gz clean_data
#
featureCounts -T 12 -p -t exon -g gene_id -a ./ref/genome.gtf -o raw_data.txt *.bam
#
rm *.bam
mkdir quick_results
mkdir ./quick_results/hisat2_log
cp *.log ./quick_results/hisat2_log
rm *.log
#
grep -v "#" raw_data.txt | awk '{for(i=1;i<=NF;i++) if ((i==1)||(6<=i)) printf("%s ",$i);print ""}'| cat | tr -s ' ' '\t' > read_counts.txt
grep -v "#" raw_data.txt | awk '{for(i=1;i<=NF;i++) if ((i==1)||(7<=i)) printf("%s ",$i);print ""}'| cat | tr -s ' ' '\t' > genes.matrix.txt
#
mv raw_data.txt.summary ./quick_results
mv raw_data.txt ./quick_results
mv read_counts.txt ./quick_results
mv genes.matrix.txt ./quick_results
rm *.bam.indel.vcf
rm *.bam.junction.bed
#end
运行脚本
run_DIY_NoRef.sh
第一步还是根据样本信息分配文件夹，随后逐个文件夹执行脚本即可
#
cat sample_list.txt | awk '{print $2}' | sort | uniq > name_list.txt
cat name_list.txt | while read line
do
mkdir ${line}_dir
cat sample_list.txt | awk -v spname=$line '{if($2 ==spname) {print $1}}' | while read line2
do
mv ${line2}* ./${line}_dir
done
done
#
cat name_list.txt | while read line
do
cp run_DIY_trinity.sh ./${line}_dir
cp run_DIY_subread.sh ./${line}_dir
cd ${line}_dir
bash run_DIY_trinity.sh
bash run_DIY_subread.sh
cd ../
done
