Bioinformatics
Blast
- 我在Linux Mint中进行数据分析，使用linuxbrew和conda安装生物软件，用brew install blast和conda install blast安装了好多次，都会报错，完全没法用，后来从Blast的FTP网站下载的预编译版本就用的好好的。看来有时还是得自己手动下载安装分析软件才行。
- 使用本地版的Blast主要分三步：安装，构建本地数据库（makeblastdb）和blastn/x/p等等
- makeblastdb 时最好带上-parse_seqids和-hash_index这两个参数
- 核苷酸比对到核苷酸数据库的常用代码：
blastn -db nt -query nt.fasta -out result.out -remote，此外还有-task参数可选，主要有三种模式：- mega-blast 比较较近源物种时使用
- blastn 比较亲缘关系较远的物种时使用
- dc-megablast 比较亲缘关系很远的物种时使用
- 其他重要的软件：diamond，blat，FAST（fastal，fastdb）
- 其他重要的参考资料：Wei Shen's Note
-
E值的公式：
- k和lambda与数据库和算法有关
- m 代表query长度
- n 代表数据库大小
- s 代表序列同源性（越大越好）
通常认为E小于10^-5就是对应S值比较可靠的结果了。
在一个数据库中E=0.001时，如果有1000条序列和querry相比对的S值比现在这个S值高，那么当E为10^-6时就只有一条序列（subject）的S值比目前这个S值更高了。
E存在的问题
- m过小时，因为S值偏小，得到的E会偏大
- 两序列的同源性高，但是有较大的Gap时，这是会有S值偏低的情况，可以通过调整Gap scores改善S值
- 有些序列的非功能区有较低的随机性，可以造成较高的同源性。
E的总结
- E适合有一定长度，且复杂度不能太低的序列
- E < 10^-5时，两序列有较高的同源性
- E < 10^-6时，两序列同源性很高，无需再次确认
Diamond
- -d
- -q
- -o
- -f
- -k
- -e
- -subject-cover
fasta 格式编辑软件
- fasta_formatter 可以用于调整fasta的格式
- bbmap的reformat.sh可以用于fasta和fastq的格式转换和多行与单行格式的互换
- idba_ud de fq2fa 可以做简单的fq到fa的格式转换，加上--merge可以实现interleaved功能，加上--filter可以实现去除包含N的功能
- bedtools getfasta可以根据gff和bed文件提取FASTA序列
vsearch软件使用
- 参数介绍：
- --iddef 指代计算序列间相似度的方式，我目前用1，仅仅考虑比对上区域的identity
- --id 指代相似度的阈值
- --masout 把cluster放到一个fasta文件中，每个cluster都包含其多序列比对和一条consensus
- --clustr-fast filename
- --cluster-smallmem 和下一个参数配合使用，可以自定义初始排序由使用者自定义
- --usersort 在cluster中，排序时很重要的，因为cluster一半用的时懒惰算法，在第一次遇见一个全新的序列时就为它开辟一个新的cluster，后面每一条序列都只和前面已经检索过的序列进行比较（我还不是很确定，所以还得再看看usearch的原理）
- --threads 指定所用的cpu数目
- --cluster string 将每个cluster分开保存，并加上由string参数指定的前缀
- --strand both 分析正向和反向序列
- 注意 ：指定输出的参数只选一个！上述加粗的参数用于设置输出格式
FAST使用
- 在将3.1GB的数据集（fastq）比对到156.6MB的database时，使用1min40s完成单核运算的步骤，当使用22个cpu核心计算剩余数据时，耗时5min20s，总物理时间7min
IGV
[to do]
- 这是一个可视化variant calling的软件
GATK
[to do]
Trimmomatic 质控软件使用
- Hiseq和Miseq系列多用Truseq3的接头
- 现在的质量标准都是Phred33，用-phred33指定，没有指定时会自动判断，但是：
the prior is phred64 - Palindrom模式的接头文件设置
- >Prefix /1
代表P5端的接头，使用Neb建库试剂盒时，为尿嘧啶U3端的碱基 - >Prefix /2
代表P7端的接头，使用Neb建库试剂盒时，是尿嘧啶U5端的碱基
- >Prefix /1
- ILLIMINACLIP：adapter.fa：最大允许去接头操作时的错配数目（2）：Palindrom去接头时的最低scores，关于score的计算方式见下文（引自Trimmomatic manual)：
Each matching base increases the alignment score by 0.6, while each mismatch reduces the alignment score by Q/10. By considering the quality of the base calls, mismatches caused by read errors have less impact. A perfect match of a 12 base sequence will score just over 7, while 25 bases are needed to score 15. As such we recommend values of between 7 - 15 as the threshold value for simple alignment mode
- LEADING：网上有些人错误理解这个参数了，这个参数设置的是序列5端最低可以接受质量，TRAILING类似，指定3端的相应参数
- CROP：指序列最大长度，超过部分从3端截掉，HEADCROP指无论如何都要从5端去掉的碱基数目
- AVGQUAL：指序列平均质量值的阈值
idba_ud
- --mink 指定最小的kmer值
- --maxk 指定最大的kmer值
- --step 指定kmer每一步增加的值
简单易用的小代码
echo your.fasta | rev | tr ATGC TACG
grep -c '^>' 可以用于计数fasta序列数目
cutadapt
[to do]
- 这是一个多才多艺的质控软件，有空再学
bioawk
[这个是神器]
- 支持sam，bed，gff，vcf和fastx格式
- 查看帮助文件非常简单 bioawk -c help
seqtk
- seq 格式转换
- comp 碱基组成信息
- sample 随机取样
- subseq 取子集
- trimfq
megahit
- --kmin
- --k-step
- -1 -2 -m -t -o --continue
- --k-list
NCBI
- 关键字得好好学学啊，例如branchiopoda[orgn]，就是检索数据的organism那个filed中包含branchiopoda关键字的record
bbmap
[这个是神器]
-
mid=97 mo=300 t=4 sort=t
这个软件可以直接在geneious中使用
FANse2
这是暨南大学张弓教授团队开发的序列比对软件，使用32bit的mpich2执行并行计算。
- Maximum read length：限制最大序列长度
- Maximum error：比对的最大错配数
- Mask genome：把参考基因组中所有的小写碱基都给屏蔽掉
- Best position：
- Min seed length：越短允许越多的error（测序错误），14比较合适
- Unidirection：关闭之后会mapping正向和反向序列
- Memory reduction：一个逻辑核心大概使用2GB当然RAM，一般情况可以不用选，选了之后大约降低20%的性能
- Export all best matches：设置在有多个等质量的比对时，是否输出多个比多结果
- Ref.Genome：顾名思义，指定参考基因组
- Dataset：待分析数据集
- Output：输出文件
- Parallel：并行核数
- split reference seg size：把待分析文件分成小块
- work folder：临时文件存储位置
demultiplex的软件
根据barcode拆分NGS data的软件有很多，可以试试seqtk_demultiplex和fastq_multx
brew
- brew cask install 与brew install的安装方式有一点点区别，前者主要暗转GUI软件，并且时预编译好了的，软件安装仅仅执行下载，解压和添加到PATH的过程，而后者会有编译和安装的过程，建议使用者先用brew cask安装，因为一般来说这样的安装在卸载时可以很干净
- brew options formula 可以查看带安装软件有什么可选参数，比如安装mrbayes时，就可以看到可以指定mpi和BEAGLE参数
- homebrew/science 不再被建议使用了，大家可以tap brewsci/bio和brewsci/science这两个formula集合
cpan
- 这是一个Perl程序的名字，让使用者可以从CPAN上下载、安装、更新以及管理在CPAN上的Perl程序
值得一看的资源
日常生活
windows与Mac共享文件
- windows访问Mac：打开我的电脑，映射网络驱动器，输入Mac的IP地址即可
- Mac访问windows：command + K，输入smb://IP，随后输入用户名和密码
linux shell
- grep技巧
- grep -w 仅仅匹配完整的单词（模式字符串左右都需要包含空格）
- grep --color=auto 高亮匹配
- -i 忽略大小写
- -e 可以指定多个匹配样式
- -f 利用文件中的模式进行匹配
- --include
- --exclude
- --exclude-from
- -B 返回匹配上面的多行
- -A 返回匹配下面的多行
- -l 返回包含匹配的文件名
-
less技巧
- 空格代表前进
- b 代表后退
- g 会到开始处
- G 去到末尾
- q退出
-
ls技巧
- -F /代表文件夹，*代表程序
- -d 屏蔽文件夹内部信息
- -r 最新的文件或文件夹放在最下面
-
tr 技巧
- tr 'a-z' 'A-Z' 小写转大写
-
uniq 技巧
- uniq -c 计数每个字符串出现的次数
-
qsub 计算机集群任务配置文件
- -I nodes=4:ppn=20,walltime=4:00:00,mem=150G
- -o fastalrun (the output dir)
- -N (the name of this process)
- $PBS_O_WORKDIR (the dir where this process was upload)
- -d (the working directory)
- -q (the priority level of this process)
-
tmux to do
-
Top 技巧（linux mint）
t查看cpu view，m查看memory view，f/F调整fileds，z：颜色，b：加粗，u/U选择带查看的用户，n显示进程数目 -
sort 技巧
-u
-t 指定间隔符号 -
nohup 服务器操作小帮手
nohup command > log 2>&1 &- 上述代码可以让你在exit登出或者是意外掉线后，远程服务器上的程序继续稳稳当当的运行，但是tmux可以做到更强大，唯一的问题就是需要服务器安装了tmux才行，然而我们大多数时候都没有办法sudo，不可以安装软件
-
zip和unzip
- zip -r file.zip file/
- unzip -t 测试文件是否完整
-
tar
- tar -xf file.tar.gz
- tar -zxvf file.tar.gz
- tar -tf file.tar.gz
-tar -ztvf file.tar.gz
f参数需要放在最后，即靠近待处理文件，t参数可以用于检视压缩文件所包含的内容，但是不解压
-
sed
- sed '2,5d' 删除第2~5行
- sed -n '5,7p' 打印出第5～7行的内容
- sed '/^$/d' 删除空白行
-
基本概念：
- 使用export的变量是环境变量，没有使用export的是自定义变量
-
scp 文件传输工具
-
split
- -a 后缀长度
- -l 行数
- file 待分解文件
- prefix 分出来的文件块的前缀
网友评论