CC BY 4.0 (除特别声明或转载文章外)
如果这篇博客帮助到你,可以请我喝一杯咖啡~
title: “生信常用的一些命令”
tags: [“生信”,”实用”,”命令”]
[TOC]
1.序列处理
1.1 基因组各种序列格式的转换:
sam 转 bam:
samtools view -@ 8 -b SAMPLE.sam > SAMPLE.bam
合并bam文件:
samtools merge -@ 30 merged.bam file1.bam file2.bam
bam 转 fastq:
samtools sort -n SAMPLE.bam -o SAMPLE_sorted.bam #排序
samtools fastq -@ 30 SAMPLE_sorted.bam > SAMPLE_sorted.fastq
fastq 转 fasta:
seqtk seq -a input.fastq > output.fasta
1.2 基因组端粒快速查询
seqtk telo -m "CCCTAA" sample.fa >> tole_info.txt
1.3 指定抽取替换基因组中的染色体
# 提取原始基因组中除了 chr30 之外的所有染色体
seqkit grep -v -p "chr30" original_genome.fa > genome_no_chr30.fa
# 把新 chr30 添加到去掉 chr30 的基因组中
cat genome_no_chr30.fa tgsgapcloser.scaff_seqs > genome_chr30_replaced.fa
2. NCBI、GSA 数据下载
2.1 iseq批量数据下载
iSeq的使用非常简单,就直接给它各大数据库的accession号就可以了,下面我介绍几个。不过,iSeq也给出了非常多的使用例子以供参考。
2.1.1 直接把ENA/SRA数据库一个项目里面的数据批量全部下载下来
iseq -i PRJNA211801
2.1.2 直接把GSA数据库一个项目里面的数据批量全部下载下来
iseq -i CRR1115842
2.1.3 直接获得单个样本的fastq文件
:
GSA数据库中只能直接获取gz压缩的fastq文件或者bam等文件,SRA/ENA数据库可以选择-q获得单个样本的fastq文件,或者-g直接下载样本的fastq.gz文件
iseq -i SRR1178105 -q
2.2 NCBI sra 测序数据下载:
下载SRA Toolkit工具:
使用命令:
prefetch SRR1482463 -O output #output替换为你想下载数据的路径
批量下载: 选中你想下载的数据,点击”Accession list”,会下载一个包含选中数据SRR号的文件(SRR_Acc_List.txt)
prefetch --option-file SRR_Acc_List.txt --max-size unlimited -O ./sra_download/
sra 转换为fastq :
fasterq-dump SRR12451710.sra -O ./fastq_output --split-files -e 8
sra 转换为fasta :
fasterq-dump --fasta SRR12451710.sra -O ./fastq_output --split-files -e 8
说明:
- -O ./fastq_output:指定输出目录;
- –split-files:将 paired-end 拆成 _1.fastq 和 _2.fastq;
- -e 8:使用 8 个线程加速转换。
3.基因组文件染色体处理
3.1 挑选指定染色体
可以使用samtools工具:
# 如果你只想保留某几条染色体,比如 Chr1 和 Chr2:
samtools faidx genome.fa Chr1 Chr2 > subset.fa
可以批量的操作:
# 先建立索引
samtools faidx genome.fa
# 列出所有序列名称
cut -f1 genome.fa.fai > all_chr.txt
# 假设要剔除 ChrM
grep -v "ChrM" all_chr.txt > keep_chr.txt
# 提取保留的染色体
samtools faidx genome.fa $(cat keep_chr.txt) > cleaned.fa
3.2 剔除指定的染色体
seqkit 支持直接剔除:
# 删除 ChrM
seqkit grep -v -p "ChrM" genome.fa > cleaned.fa
多个染色体一起剔除:
seqkit grep -v -r -p "Chr2|Chr1|Chr3" genome.fa > cleaned.fa
3.3 截取基因组中某条染色体
如果有一条序列 genome.fa,比如有 Chr1,你想在 1,000,000 bp 和 2,000,000 bp 处分割:
# 先建立索引
samtools faidx genome.fa
# 截取 1~1,000,000
samtools faidx genome.fa Chr1:1-1000000 > Chr1_part1.fa
# 截取 1,000,001~2,000,000
samtools faidx genome.fa Chr1:1000001-2000000 > Chr1_part2.fa
# 截取 2,000,001 到末尾
samtools faidx genome.fa Chr1:2000001- > Chr1_part3.fa