宏基因组:MAG丰度计算(个人想法)

步骤:将去冗余的MAGs的所有contig,构建bin_contig关系后,合并为一个fasta文件,以这个文件为reference,构建bwa-mem2的index,然后将所有样品的clean_reads映射map到这个reference上,根据比对情况求取所有contigs的TPM,最后将同一bin的contigs的TPM进行求和,得到MAGs的TPM表格

# 构建bin_contig关系
parse_stb.py --reverse -f bin/* -o bin/goat.stb
# 合并所有contigs
sed 's/>/\n>/g' bin/*fa | sed '/^$/d' > bin/797MAGs.fa
# 构建index
/bwa-mem2/bwa-mem2 index bin/797MAGs.fa

for i in *_2.fastq
do
	num=${i%%_2.fastq}
	/bwa-mem2/bwa-mem2 mem -t 180 \
	bin/797MAGs.fa ${num}_1.fastq ${num}_2.fastq > bin/${num}_aligen.sam
	samtools view -Sb -@180 bin/${num}_aligen.sam > bin/${num}_aligen.bam
	samtools sort -@180 bin/${num}_aligen.bam -o bin/${num}_sort.bam
	samtools index -@180 bin/${num}_sort.bam
	samtools depth bin/${num}_sort.bam > bin/${num}_depth.txt
	samtools idxstats -@180 bin/${num}_sort.bam > bin/${num}_aligen_result.txt
	sed -i "1 i GeneID\t\length\tmapped_read\tunmapped_read" bin/${num}_aligen_result.txt
# 计算TPM
	/my_script/TPM_count.py bin/ ${num}_aligen_result.txt ${num}_TPM.txt
# 计算coverage和depth
	/my_script/coverage_depth_filter.py --ref_fa bin/797MAGs.fa --depth_file bin/${num}_depth.txt --output_file bin/${num}_cunzai.txt --coverage 0 --average_depth 0

	rm bin/${num}_aligen.sam
	rm bin/${num}_aligen.bam
	rm bin/${num}_sort.bam
	rm bin/*bam*
	rm bin/${num}_depth.txt
done

# 把各样品的TPM合并为一个表
/my_script/TPM_together.py --work_path . --file_maker _TPM.txt --output_name 797_MAGs_tpm.txt

# 根据bin把这个表求和
/my_script/TPM_together_by_bin_contig.py --input_file 797_MAGs_tpm.txt --input_bin_contig 797_goat.stb --output_file 797_MAGs_tpm_by_bin.txt
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值