相关知识幕布整理
https://mubu.com/doc/2tQ_CfEmMg
基础20题
题目链接：http://www.bio-info-trainee.com/2900.html
#级联创建文件夹
mkdir -p 1/2/3/4/5/6/7/8/9
#添加文件，输入内容，ctrl+d退出
cat > me.txt
#查看文件内容
cat me.txt
#删除以上文件以及文件夹
rm -r 1
#批量创建文件夹
mkdir -p folder{1..5}/folder{1..5}
#不同文件夹下都创建同一文件
cat > me.txt
echo folder{1..5}/folder{1..5} |xargs -n1 cp me.txt
#删除建立的文件以及文件夹
rm -r folder{1..5}
#下载文件
wget http://www.biotrainee.com/jmzeng/igv/test.bed
##查看文件行数
wc -l test.bed
##查看关键词在第几行
grep -n "H3K4me3" test.bed
#下载压缩文件并解压
wget http://www.biotrainee.com/jmzeng/rmDuplicate.zip && unzip *.zip
#查看文件夹结构
tree rmDuplicate
#打开bam文件
samtools view -H tmp.sorted.bam
#计算bam文件第二列0和16数字个数
samtools view tmp.sorted.bam | cut -f2 | sort -n | uniq -c
#同样解压，统计以 >> 开头有多少行
cat fastqc_data.txt | grep '>>' | wc -l
#在hg38.tss文件中找基因，如OTOF基因
cat hg38.tss | grep 'NM_001287489'
##结果为 NM_001287489 chr2 26556698 26560698 1
#每条染色体的基因个数——染色体名出现次数
cut -f2 hg38.tss | sort -rn | uniq -c
#统计NM、NR开头的序列
##NM ： 转录产物序列，成熟mRNA转录本序列
##NR ： 非编码的转录子序列，包括结构RNAs，假基因转子
cat hg38.tss | grep 'NM' | wc -l
cat hg38.tss | grep 'NR' | wc -l
数据格式习题
fasta和fastq练习
fasta格式是一种基于文本的，用于表示核苷酸序列或氨基酸序列的格式。第一行是由'>'开头的任意文字说明，用于序列标记，第二行开始为序列本身。
fastq格式是基于文本的，保存核酸序列和其测序质量信息的标准格式，四行为一个基本单元，第一行记录序列标识以及相关的描述信息，以'@'开头，单个序列标识具有唯一性，第二行为碱基序列，第三行以'+'开头，第四行是质量信息，长度与第二行相对应。
题目链接：http://www.bio-info-trainee.com/3575.html
1、统计reads_1.fq 文件中共有多少条序列信息
cat reads_1.fq |wc -l
# 除以4就是序列数
2、输出所有的reads_1.fq文件中的标识符(即以@开头的那一行)
cat reads_1.fq | paste - - - - | cut -f1
# paste - - - - paste命令用于合并文件的列，这里将每四行变成四列
3、输出reads_1.fq文件中的 所有序列信息(即每个序列的第二行)
cat reads_1.fq | paste - - - - | cut -f2
# 用awk提取 awk '{if(NR%4==2)print}' reads_1.fq 思路同
4、输出以‘+’及其后面的描述信息(即每个序列的第三行)
cat reads_1.fq | paste - - - - | cut -f3
5、输出质量值信息(即每个序列的第四行)
cat reads_1.fq | paste - - - - | cut -f4
# 发现awk要这样才可以：awk '{if(NR%4==0)print}' reads_1.fq
6、计算reads_1.fq 文件含有N碱基的reads个数
awk '{if(NR%4==2)print}' reads_1.fq |grep -c N
# grep -c 只输出匹配行的计数
7、统计文件中reads_1.fq文件里面的序列的碱基总数
awk '{if(NR%4==2)print}' reads_1.fq |grep -o [ATGCN]|wc
# grep -o 只输出文件中匹配到的部分
8、计算reads_1.fq 所有的reads中N碱基的总数
awk '{if(NR%4==2)print}' reads_1.fq | grep -o N |wc
9、统计reads_1.fq 中测序碱基质量值恰好为Q20的个数
10、统计reads_1.fq 中测序碱基质量值恰好为Q30的个数
awk '{if(NR%4==0)print}' reads_1.fq | grep -o 5 |wc
awk '{if(NR%4==0)print}' reads_1.fq | grep -o ? |wc
11、统计reads_1.fq 中所有序列的第一位碱基的ATCGNatcg分布情况
awk '{if(NR%4==2)print}' reads_1.fq | cut -c 1 |sort|uniq -c
# -c 1 切分排序去重
12、将reads_1.fq 转为reads_1.fa文件(即将fastq转化为fasta)
cat reads_1.fq | paste - - - - | cut -f1,2|tr '\t' '\n'|tr '@' '>' > reads_1.fa
# 开头符号变化
13、 统计上述reads_1.fa文件中共有多少条序列
cat reads_1.fa |wc -l
# 除以2
14、计算reads_1.fa文件中总的碱基序列的GC数量
cat reads_1.fa |grep -o [GC] |wc -l
15、删除 reads_1.fa文件中的每条序列的N碱基
cat reads_1.fa |tr -d "N"
# tr -d 删除
16、删除 reads_1.fa文件中的含有N碱基的序列
cat reads_1.fa |paste - -|grep -v N | tr '\t' '\n'
17、删除 reads_1.fa文件中的短于65bp的序列
cat reads_1.fa |paste - -|awk '{if (length($2)>65) print}'
18、删除 reads_1.fa文件每条序列的前后五个碱基
head reads_1.fa|paste - - | cut -f2|cut -c 5-
19、删除 reads_1.fa文件中的长于125bp的序列
cat reads_1.fa | paste - -|awk '{if (length($2)<125) print}'| wc
20、查看reads_1.fq 中每条序列的第一位碱基的质量值的平均值
awk '{if(NR%4==0)print}' reads_1.fq | cut -c 1 |
sam和bam练习
SAM是一种序列比对格式标准，由sanger制定，是以TAB为分隔符的文本格式。主要应用于测序序列mapping到基因组上的结果表示，也可以表示多重比对的结果。SAM分为两部分，注释信息(header section)和比对结果部分(alignment section)。
BAM是SAM的二进制格式，占用储存空间更小，运算更快。
题目链接：http://www.bio-info-trainee.com/3578.html
# less -S tmp.sam 查看
cat tmp.sam |grep -v "^@"|wc
# 共有多少种比对类型
cat tmp.sam |grep -v "^@"|cut -f 2|sort|uniq -c
# 筛选比对失败的reads，第六列的*代表失败，特征为N数量多
cat tmp.sam | grep -v '^@' |awk '{if($6=="*")print}'|wc
# 比对失败的reads区分成单端失败和双端失败情况，并且拿到序列ID
# 单端失败
cat tmp.sam | grep -v '^@' |awk '{if($6=="*")print $1}'|sort|uniq -c|grep -w 1
# 双端失败
cat tmp.sam | grep -v '^@' |awk '{if($6=="*")print $1}'|sort|uniq -c|grep -w 2
# 筛选出比对质量值大于30的情况，看第五列
cat tmp.sam | grep -v '^@' |awk '{if($5>30)print}'|wc
# 筛选比对成功，但不完全匹配的序列
cat tmp.sam | grep -v '^@' |awk '{if($6!="*")print$6}'|grep "[IDNSHPX]"|wc
# 筛选出inset size长度大于1250bp的 pair-end reads
cat tmp.sam | grep -v '^@' |awk '{if($7>1250)print}'|less -S
# 统计参考基因组上面各条染色体的成功比对reads数量
cat tmp.sam | grep -v '^@' |cut -f 3|sort -u
# 筛选出原始fq序列里面有N，但是比对的时候却是完全匹配的情况，第十列匹配N
cat tmp.sam | grep -v '^@'|awk '{if($10~N)print}'|awk '{if($6!~"[IDNSHP]")print}'|awk '{if($6!="*")print}'|less -S
# sam文件里面的头文件行数
cat tmp.sam |grep "^@"|wc
# sam文件里每一行的tags个数一样吗?不一样
# sam文件里每一行的tags个数分别是多少个？
运用awk数组
# sam文件里记录的参考基因组染色体长度分别是？
head tmp.sam | grep 'LN'
# 找到比对情况有insertion情况的；找到比对情况有deletion情况的
cat tmp.sam | grep -v '^@'|awk '{if($6~I)print}'|less -S
cat tmp.sam | grep -v '^@'|awk '{if($6~D)print}'|less -S
# 取出位于参考基因组某区域的比对记录，比如 5013到50130 区域
cat tmp.sam | grep -v '^@'|awk '{if($4>5013 && $4 <50130)print}'|less -S
# 把sam文件按照染色体以及起始坐标排序
cat tmp.sam | grep -v '^@'|awk '{print $4}'|sort -n
# 找到 102M3D11M 的比对情况，计算其reads片段长度
grep 102M3D11M tmp.sam|cut -f 10|wc
VCF练习
VCF，Variant Call Format 存储变异位点的标准格式，可以用来表示单核苷酸多态性SNP、插入缺失、结构变异、拷贝数量变异。VCF使用UTF-8编码，以##开头的注释信息和具体突变信息。
题目链接：http://www.bio-info-trainee.com/3577.html
# 把突变记录的vcf文件区分成 INDEL和SNP条目
cat tmp.vcf | grep -v '^#'|less -S|awk '{if (length($4)==1 && length($5)==1) print}'
# 统计INDEL和SNP条目的各自的平均测序深度
grep -v '^#' tmp.vcf |cut -f 8|head -10|grep 'DP='
# 把INDEL条目再区分成insertion和deletion情况
cat tmp.vcf | grep -v '^#'|less -S|awk '{if (length($4)<length($5)) print}'
cat tmp.vcf | grep -v '^#'|less -S|awk '{if (length($4)>length($5)) print}'
# 统计SNP条目的突变组合分布频率
cat tmp.vcf |grep -v '^#'|awk '{if (length($4)==1 && length($5)==1) print}'|cut -f4,5|sort|uniq -c
# 找到基因型不是 1/1 的条目，个数
cat tmp.vcf |grep -v '^#'|awk '{ print $10}'|grep -v '^1/1'
# 筛选测序深度大于20的条目
cat tmp.vcf |grep -v '^#' |awk '{if($11>20);print}'
# 筛选变异位点质量值大于30的条目
cat tmp.vcf |awk '{if ($6>30) print}'|grep -v '^#'|less -S|head
# 组合筛选变异位点质量值大于30并且深度大于20的条目
cat tmp.vcf |awk '{if ($6>30 && $11>20) print}'
好久之前做的练习，现在Linux算是跨过门槛了，但是不熟练，Linux、R、Python都是我们处理数据的工具，需要的时候能用它们处理问题，是我目前对自己编程方面的要求啦~
更多学习资源：
生信技能树公益视频合辑
生信技能树简书账号
生信工程师入门最佳指南
生信技能树全球公益巡讲
招学徒
...
你的宣传能让数以万计的初学者找到他们的家，技能树平台一定不会辜负每一个热爱学习和分享的同道中人 😎
