生信数据格式转换汇总（一）

zhendehaibucuo

已于 2024-07-08 09:45:53 修改

阅读量1k

点赞数 14

文章标签： linux

于 2024-07-05 16:47:18 首次发布

本文链接：https://blog.csdn.net/zhendehaibucuo/article/details/140211006

版权

fasta格式转换

由于一些软件对输入的fasta文件格式具有要求有时需要将其从一行转化为多行，有事需要将其从多行转化成一行。这时候我们就需要用到FASTX-Toolkit工具将其进行转化。

FASTX-Toolkit下载

wget http://hannonlab.cshl.edu/fastx_toolkit/fastx_toolkit_0.0.13_binaries_Linux_2.6_amd64.tar.bz2
tar xjvf fastx_toolkit_0.0.13_binaries_Linux_2.6_amd64.tar.bz2

主要参数

(haphic) [liyanchun@login1 00.data]$ fasta_formatter -h
usage: fasta_formatter [-h] [-i INFILE] [-o OUTFILE] [-w N] [-t] [-e]
Part of FASTX Toolkit 0.0.14 by assafgordon@gmail.com

   [-h]         = This helpful help screen.
   [-i INFILE]  = FASTA/Q input file. default is STDIN.
   [-o OUTFILE] = FASTA/Q output file. default is STDOUT.
   [-w N]       = max. sequence line width for output FASTA file.
                  When ZERO (the default), sequence lines will NOT be wrapped -
                  all nucleotides of each sequences will appear on a single 
                  line (good for scripting).
   [-t]         = Output tabulated format (instead of FASTA format).
                  Sequence-Identifiers will be on first column,
                  Nucleotides will appear on second column (as single line).
   [-e]         = Output empty sequences (default is to discard them).
                  Empty sequences are ones who have only a sequence identifier,
                  but not actual nucleotides.

主要用到的参数就是-i、-o、-w。前俩不用过多介绍即输入输出。-w为主要参数将其设置为-w 0 即将多行转化为一行；若将其设置为-w 60即一行转化为多行，且每行有60个碱基。

示例

这是转化前用命令 less -S asm.fa 打开的asm.fa文件格式

fasta_formatter -i asm.fa -o sub.fa -w 0

使用该命令转化后用 less -S sub.fa 打开sub.fa文件

fasta_formatter -i sub.fa -o sub20.fa -w 20

再使用该命令将其转化为每个20个碱基的表达形式

fastq转化为fasta

cat gene.fastq |paste - - - - |cut -f 1,2|sed 's/^@/>/'|tr "\t" "\n"|awk '{print $1}'>gene.fa

大家在使用时候只需要将gene.fastq换为自己的文件即可，这个功能其实好多工具都可以实现，大家没必要纠结，包括后续要介绍到的seqtk也可以实现。seqtk的功能非常强大，我们后续也会用到，这里就不过多介绍感兴趣的同学可以自行了解。

冲bam文件中提取到fastq文件

samtools sort -@ 20 -n your.bam  | samtools fastq -@ 20  > output.fastq

这里我们选择使用samtools工具进行提出，-@ 20 即线程数设置为20；-n your.bam 即为你自己的bam文件；output.fastq即为输出的fastq序列名字。

SRA文件转化为fastq文件

有时从NCBI上下载下载下来的文件并不是bam文件而是sra文件，这时候我们就需要用到SRA_Toolkit工具就将其转化为fastq文件。

fastq-dump --split-files SRR9331727

注意这时候就有一个小小滴坑点来喽，当你的文件时Hi-C等双端测序序列时候使用此命令，当你时带三代单端测序数据时直接使用

fastq-dump  SRR9331727

另外，补充一下从NCBI上批量下载sra序列和批量转化的方法

#下载
cat my_Acc_List.txt | while read id;do (prefetch -X 100G $id );done
#转化
ls -d SRR* | while read id;do ( fastq-dump  $id  );done

创建my_Acc_List.txt并将自己需要下载序列的序列号输入即可。

gfa转化为fa

grep "^S" p_utg.gfa | awk '{print ">"$2"\n"$3}' > p_utg.gfa.fa

sam转化为bam

samtools view -bS sd_0016.sam > sd_0016.bam

提取特定名称的序列

#!/usr/bin/python
import sys
import gzip


def get_chr_len(in_fasta, out_file, is_chr_only):
	if is_chr_only != "T" and is_chr_only != "F":
		print("Error argument")
		exit(0)
	dict_len = {}
	if in_fasta[-3:].lower() == '.gz':
		f_in = gzip.open(in_fasta, 'rt')
	else:
		f_in = open(in_fasta, 'r')
	chrn = ''
	seq = ''
	for line in f_in:
		if line[0] == '>':
			if seq != '':
				dict_len[chrn] = len(seq)
			chrn = line.strip().split()[0][1:]
			seq = ''
		else:
			seq += line.strip()
	dict_len[chrn] = len(seq)
	f_in.close()
	
	with open(out_file, 'w') as f_out:
		chr_list = sorted(dict_len.keys())
		for chrn in chr_list:
			if is_chr_only == "T" and chrn[:3].lower() != 'chr':
				continue
			else:
				f_out.write(chrn+"\t"+str(dict_len[chrn])+"\n")


if __name__ == "__main__":
	if len(sys.argv) < 4:
		print("Notice: script for calculating length of chromosomes in fasta file")
		print("Usage: python "+sys.argv[0]+" <fasta_file> <output_file> <T/F chr only>")
	else:
		f_fasta = sys.argv[1]
		f_out = sys.argv[2]
		is_chr_only = sys.argv[3]
		get_chr_len(f_fasta, f_out, is_chr_only)

这里补充一个提取序列长度的脚本，并将其命名为get_chr_len.py。

python get_chr_len.py scaffolds_FINAL.fasta 1_length F#(1)改用该脚本获得每条序列的长度
sort -k 2 -r -n 1_length > 1_length.sort#(2)对长度进行排序
awk 'NR <= 24 {print $1}' 1_length.sort > list#（3）提取最长的24条序列的名称
awk -F'[ >]' 'NR==FNR {wanted[$1]; next} /^>/ {found=($2 in wanted)} found' list scaffolds_FINAL.fasta > hap_hifi_yash.fasta#(4)将其提取出来

当然你可能并不想提取最长的几条序列，这并不影响，只需要将你需要提取的序列命令输入list直接使用第四步提取即可。

除此之外还可以使用seqtk进行提取。

seqtk subseq scaffolds.fa name.txt >sub.fa

注意！！！！！！！

上述提到的所有保存名字的文件均需以以下形式保存（一行一个！！！！）

还是老规矩，如果本文中有哪些命令看不懂的找老朋友gpt，当然如果就想让我解答可以留言或者给我发私聊，欢迎大家提问。当然如果大家有其他想要了解的数据格式转化的方法，或者有想要分享的方法欢迎大家在评论区留言哦~

参考：

FASTX-Toolkit - Command Line Usage

[Linux|生信]project4_01:批量下载sra文件并转化为fastq文件_sra批量下载转化fastq-CSDN博客