宏基因组数据 数据预处理

下载数据

# 需提供SRR号
# 下载sra文件使用prefetch SRR文件 -p -o 自定义名称（这个自定义数据为之后的$1 $2数据）默认即可
prefetch SRR1262938 SRR1263009 -p # -p显示进度

拆分双端数据

# 输入文件为上面的.sra文件 
# -p 显示进度 -e 选择cpu线程
fasterq-dump -3 SRR1262938 -p -e 100 -O out_dir # out_dir为输出文件夹

注：这是shotgun的处理过程，HiC的处理过程和上面相同，最终结果如下

处理原始read

# 预处理 三步应该可以放一块吧，别人论文这样写的我就这样弄了
# step1 ./bbduk.sh  in1=() in2=() out1=() out2=() ref=adapters.fa ktrim=r k=23 mink=11 hdist=1 minlen=50 tpe tbo
./bbduk.sh  in1=SRR1262938_1.fastq in2=SRR1262938_2.fastq out1=out_SRR1262938_1.fastq out2=out_SRR1262938_2.fastq ref=/home/yanziming/vicent/bbmap/resources/adapters.fa ktrim=r k=23 mink=11 hdist=1 minlen=50 tpe tbo
# step2 ./bbduk.sh  in1=() in2=() out1=() out2=() trimq=10 qtrim=r ftm=5 minlen=50
./bbduk.sh  in1=out_SRR1262938_1.fastq in2=out_SRR1262938_2.fastq out1=out1_SRR1262938_1.fastq out2=out2_SRR1262938_2.fastq trimq=10 qtrim=r ftm=5 minlen=50
# step3  ./bbduk.sh  in1=() in2=() out1=() out2=() ftl=10
./bbduk.sh in1=out1_SRR1262938_1.fastq in2=/home/yanziming/vicent/data_set/synthec_metagenomic_yeast/shotgun/SRR1262938/out2_SRR1262938_2.fastq out1=finre_SRR1262938_1.fastq out2=finre_SRR1262938_2.fastq ftl=10

组装shotgun（两个软件）

meaghit

# $shotgun_file1 $shotgun_file2为上面剪枝处理后文件
megahit -1 $shotgun_file1 -2 $shotgun_file2 -o 'megahit_out' --min-contig-len 1000 --k-min 21 --k-max 141 --k-step 12 --merge-level 20,0.95
#下图为输出文件夹内容 重点要用的文件是 final.contigs.fa

spades.py

# $shotgun_file1 $shotgun_file2为上面剪枝处理后文件 该方法不支持定义最小contigs长度
spades.py --meta -1 $shotgun_file1 -2 $shotgun_file2 -o 'spades_out' -t 66
# 下图为输出文件夹 重点要用的文件是 contigs.fasta

生成Hi-C图的第一个文件fa文件

bwa比对

# 建立索引 $fa就是上面提到的fa文件
bwa index $fa
# bwa比对 $hic_file1 $hic_file2 是上面处理得到的两个Hi-C文件 输出文件为MAP.sam -t 表示线程数
bwa mem -5SP -t 100 $fa $hic_file1 $hic_file2 > 'MAP.sam'
# MAP.sam为上面的输出文件 路径可以自定义
samtools view -F 0x904 -bS 'MAP.sam' > 'MAP_unsorted.bam'
# samtools排序
samtools sort -n 'MAP_unsorted.bam' -o 'MAP_sorted.bam' -@ 66

生成Hi-C图的第二个文件MAP_sorted.bam

标签生成->blastn&处理文件

# 比对 其他参数默认，-num_threads cpu核数 -qurey 上面组装得到的fa文件 -out 输出的标签文件
/media/ubuntu/conda/vicent/blast/ncbi-blast-2.2.28+/bin/blastn -query final.contigs.fa -db /media/ubuntu/conda/vicent/blast/ncbi-blast-2.2.28+/db/nt -perc_identity 95 -evalue 1e-30 -word_size 50 -num_threads 40 -outfmt "6 qacc evalue pident stitle" -out ./m_yeast.txt
# 这一步需要在九楼工作站进行

# 该脚本针对合成宏基因组，其他的数据集需要更改对应标签，
import csv
import re
from tqdm import tqdm
# for i in tqdm(range(n), desc="my_sort"):
f = open("/home/yanziming/vicent/data_set/synthetic_metagenomic_yeast/shotgun/m_yeast.txt", 'r')
content = []
lines = f.readlines()
for line in tqdm(lines, desc="处理初始数据"):
    line = line.replace('\n', '')
    line = line.split('\t')
    line[3] = ' '.join(line[3].split()[1:3])
    content.append([line[0], line[3]])
# print(content[0])
label_list = ['Saccharomyces cerevisiae', 'Saccharomyces paradoxus',
              'Saccharomyces mikatae', 'Saccharomyces kudriavzevii',
              'Saccharomyces bayanus', 'Naumovozyma castellii',
              'Lachancea waltii', 'Lachancea kluyveri',
              'Kluyveromyces lactis', 'Kluyveromyces wickerhamii',
              'Ashbya gossypii', 'Scheffersomyces stipitis',
              'Pichia pastoris']

node_name = ''
# 提取每个node标签
node_label = []
for row in tqdm(content, desc='过滤每个contigs标签'):
    if row[1] in label_list and node_name != row[0]:
        node_label.append(row)
        node_name = row[0]
# node_label存储对应的节点名和标签
for i in range(10):
    print(node_label[i])
print(len(node_label))
# 将 node_label输出为csv文件作为生成Hi-C图的TAX参数
path = '/home/yanziming/vicent/data_set/synthetic_metagenomic_yeast/process_data/label.csv'
with open(path, 'w', newline='') as csvfile:
    # newline=''参数是必需的，以确保在Windows操作系统上写入CSV文件时不会出现额外的空行
    writer = csv.writer(csvfile)
    writer.writerows(node_label)

生成Hi-C图的第三个文件label.csv

覆盖率生成

注：使用bbmap.sh文件进行比对，用shotgun read比对组装的contigs

bbmap.sh in1=SG1.fastq.gz in2=SG2.fastq.gz ref=final.contigs.fa out=SG_map.sam bamscript=bs.sh; sh bs.sh# 注意这一行是bbmap.sh的参数
jgi_summarize_bam_contig_depths --outputDepth coverage.txt --pairedContigs pair.txt SG_map_sorted.bam

生成Hi-C图的第四个文件coverage.txt

Hi-C图生成

利用上面四个文件，输入Hi-C图生成脚本中，得到Hi-C图

# 运行脚本
python hiczin_contact.py
# 将npz转为txt文件
python npz_to_txt.py

kmer特征生成

脚本在/home/yanziming/vicent/my_research_contents/my_code/node_feature/calculate_kmer.py中，用于生成不同kmer的特征向量

python calculate_kmer_multi_thread.py