組裝策略
適用于設(shè)計(jì)多樣本多物種的組裝。例如100個(gè)樣本，10個(gè)物種混聊。這里如果想直接完成10個(gè)de nove組裝嘴办，需要將所有樣本數(shù)據(jù)放到一起后瞬场，通過樣本信息表聲明每個(gè)樣本的物種名，再分別組裝涧郊。因此涉及通過awk操作樣本信息表贯被。另外，多物種的分析中涉及到直系同源基因妆艘。此時(shí)通過unigene進(jìn)行連接是最佳的彤灶。但在transdecoder預(yù)測(cè)ORF時(shí)，有時(shí)同一個(gè)轉(zhuǎn)錄本預(yù)測(cè)得到兩個(gè)不同的ORF批旺。所以最終是以預(yù)測(cè)得到的ORF為單元進(jìn)行定量和分析幌陕。響應(yīng)的CDS和PEP用以計(jì)算直系同源基因以及后續(xù)進(jìn)化分析等。因此汽煮，這個(gè)過程中搏熄，從Trinity.fasta得到unigene，再以u(píng)nigene為參考基因組暇赤，預(yù)測(cè)ORF得到gff3和gtf文件心例，進(jìn)行有參組裝。這種方式用來規(guī)避用Trinity.fasta為參考時(shí)翎卓，涉及到同一個(gè)基因得到不同注釋契邀，以及后續(xù)提取最長(zhǎng)轉(zhuǎn)錄本時(shí)遇到的麻煩等。
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
按照第二列進(jìn)行分配
aaa="SP"
cat sample_list.txt |
awk -v bbb=$aaa '{if($2 ==bbb) {print $1,$2}}'
2023NJCYX303-SZ SP
2023NJCYX302-SZ SP
事實(shí)上失暴，我們想要將第一列的樣本按文件夾分配
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按照第二列的名稱進(jìn)行循環(huán)
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
這樣所有的樣本被歸到響應(yīng)的文件夾里
Trinity組裝（單個(gè)）
基本代碼如下：
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
會(huì)報(bào)錯(cuò)坯门，反正經(jīng)常這樣報(bào)錯(cuò)：
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.
類似的報(bào)錯(cuò)一般是先手動(dòng)執(zhí)行報(bào)錯(cuò)的這一步，獲得left.fa和right.fa逗扒，后面再運(yùn)行就對(duì)了古戴。有時(shí)候會(huì)順手直接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
修改名稱
由于前面根據(jù)物種信息表矩肩，將物種名信息放在文件夾里∠帜眨現(xiàn)提取文件夾的物種名信息，加入unigene的文件名文件內(nèi)的序列中。
注意叉袍，“物種名|Trinity名”這種格式雖然很不錯(cuò)始锚，但“|”在建立索引時(shí)會(huì)報(bào)錯(cuò)，所以舍棄了喳逛。
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
預(yù)測(cè)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
我們首先有一個(gè)自動(dòng)組裝trinity的腳本典蝌，這個(gè)腳本僅完成到ORF注釋曙砂，因?yàn)椴淮_定每個(gè)流程中需要哪些步驟（未必每次都需要定量或者功能注釋）
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
# 并非重復(fù)，第一遍經(jīng)常運(yùn)行失敗
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
有參組裝腳本骏掀，實(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
運(yùn)行腳本
run_DIY_NoRef.sh
第一步還是根據(jù)樣本信息分配文件夾鸠澈，隨后逐個(gè)文件夾執(zhí)行腳本即可
#
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
