【原】【Linux 笔记】Linux 基本操作 - 04. 利用SRA数据实践文本处理

hucy_Bioinfo 2021-09-04

展开全文

笔记接上篇【Linux 笔记】Linux 基本操作 - 03. shell脚本编程。笔记大部分源于生信技能树的B站视频教程【生信技能树】生信人应该这样学linux(更新至第14集)，如有需要，可去欣赏原汁原味的视频讲解。

10. 文本处理实践

下载SRA数据库中SRP115453 (PRJNA398328) 研究的测序数据信息文件 SraRunTable.txt 和测序结果ID文件SRR_Acc_List.txt，利用FileZilla或者WinSCP将其上传到服务器，以此作为此次文本操作实践的测试数据。数据库的介绍详见：【数据库】SRA数据库介绍及数据下载

测试数据文件下载地址：
https://www.ncbi.nlm./Traces/study/?acc=SRP115453&o=acc_s%3Aa

文章：Liu H, Murphy CJ, Karreth FA, Emdal KB, White FM, Elemento O, Toker A, Wulf GM, Cantley LC. Identifying and Targeting Sporadic Oncogenic Genetic Aberrations in Mouse Models of Triple-Negative Breast Cancer. Cancer Discov. 2018 Mar;8(3):354-369. doi: 10.1158/2159-8290.CD-17-0679.
文章下载地址：http://cancerdiscovery./content/candisc/early/2018/02/14/2159-8290.CD-17-0679.full.pdf
数据地址：
https://www.ncbi.nlm./sra/?term=SRP115453
https://www./ena/browser/view/PRJNA398328?show=reads

基本操作

grep RNA SraRunTable.txt | wc  #统计RNA-seq的数据个数
head SRR_Acc_List.txt > id     #数据量较大，取前10个作为测试
cat id

# 将默认名字修改为简单易区分的名称
mv SraRunTable.txt info   
mv SRR_Acc_List.txt allid
head * | less -NS #q；查看当前路径下所有文件的前10行

grep 筛选数据（行）

grep -f id info | wc # -f:在文件中进行匹配; 
# grep的参数
# -w：完全匹配(word)
# -c：计数，grep -fc id info
# -v：返回未匹配到的部分
grep -v '^#' tmp #匹配不以#开题的行

which grep  # 所在目录 /bin/grep
type grep  # grep 是 `grep --color=auto' 的别名

cut 和 awk 操作数据（列）

less -SN info # -N显示行号；-S:不换行显示

# 查看文件的表头，按列显示，并标注行号，方便后续的操作
head -1 info | tr ' ' '\n' | cat -n  #tr: 将以空格分隔的列，修改为按换行符分隔。cat -n: 显示行号
## 1  Run,age,Assay
## 2  Type,AvgSpotLen,Bases,BioProject,BioSample,BioSampleModel,Bytes,Center
## 3  Name,Consent,DATASTORE
## 4  filetype,DATASTORE
## 5  provider,DATASTORE
## 6  region,Experiment,Genotype,Instrument,Library
## 7  Name,LibraryLayout,LibrarySelection,LibrarySource,mouse_number,Organism,Platform,primary_tumor_index,ReleaseDate,Sample
## 8  Name,sequencing_type,sex,SRA
## 9  Study,Strain,Tissue,Treatment,tumor_index,replicate_number

# 数据的存储，比想象中更复杂一些，接下来现将所有列的分隔符，由空格改成逗号分隔
head -1 info | tr ' ' ',' |  tr ',' '\n' | cat -n
## 1  Run
## 2  age
## 3  Assay
## 4  Type
## 5  AvgSpotLen
## 6  Bases
## 7  BioProject
## 8  BioSample
## 9  BioSampleModel
## 10  Bytes
## 11  Center
## 12  Name
## ……
head -2 info | tail -1 | tr ' ' ',' |  tr ',' '\n' | cat -n

head -4 info | tr ' ' ',' |  tr ',' ' ' | cut -d" " -f1,4,7,8,27,43-47  #-d: 指定分隔符
##Run Type BioProject BioSample LibrarySource replicate_number
##SRR5933808 WXS PRJNA398328 SAMN07504110 Trp53+/+;Brca1+/+; 2017-11-11T00:00:00Z HL-WES-34 WEX female SRP115453
##SRR5933809 WXS PRJNA398328 SAMN07504109 Trp53flox/flox;Brca1+/+; 2017-11-11T00:00:00Z HL-WES-22 WEX female SRP115453
##SRR5933810 RNA-Seq PRJNA398328 SAMN07504108 Trp53flox/flox;Brca1flox/flox; HL145 RNAseq female SRP115453 C57BL/6

# 可见上述第27列，是由分号分隔的元素，取分号前的一列
head -4 info | tr ' ' ',' |  tr ',' ' ' | cut -d" " -f27 | cut -d";" -f 1
##LibrarySource
##Trp53+/+
##Trp53flox/flox
##Trp53flox/flox

head -4 info | tr ' ' ',' |  tr ',' ' ' | cut -d" " -f27 | tr ";"  "\t" #替换分隔符
head -4 info | tr ' ' ',' |  tr ',' ' ' | cut -d" " -f27 | sed 's/;/\t/g' #替换分隔符
##LibrarySource
##Trp53+/+        Brca1+/+
##Trp53flox/flox  Brca1+/+
##Trp53flox/flox  Brca1flox/flox
##Trp53+/+        Brca1+/+

## 查看该研究中的实验类型，排序并统计频数
cat info | tr ' ' ',' |  tr ',' ' ' | cut -d" " -f4 | sort | uniq -c
## 88 RNA-Seq
## 94 WXS
## 可见在182个实验中，RNA-Seq数据共88个，WXS数据共94个。

# 对染色体进行排序
# ln -s /mnt/hd3/shuju/hucy/humandb/hg19_refGene.txt /mnt/hd3/shuju/hucy/SRA_data/hg19_refGene.txt
cut -f3 hg19_refGene.txt | sort | uniq -c
## 7053 chr1
## 3465 chr10
## 4125 chr11
## ……
## 3712 chr17
## 63 chr17_ctg5_hap1
## 1 chr17_gl000205_random
## 1175 chr18
## 4278 chr19
## 21 chr19_gl000209_random
## 9 chr1_gl000191_random

# 按染色体大小号顺序排序
cut -f3 hg19_refGene.txt | sort -V | uniq -c  #-V, --version-sort：在文本内进行自然版本排序
## 7053 chr1
## 9 chr1_gl000191_random
## 1 chr1_gl000192_random
## 4938 chr2
## 4256 chr3
## 2639 chr4
## 14 chr4_ctg9_hap1
## 1 chr4_gl000193_random
## 4 chr4_gl000194_random
## 3264 chr5
## 3628 chr6
## 256 chr6_apd_hap1
## 3780 chr7
## 9 chr7_gl000195_random
## 2852 chr8
## 2900 chr9
## 3465 chr10
## 4125 chr11

cut -f3 hg19_refGene.txt | sort | uniq -c | sort -k1,1nr | head  #按染色的频数排序
cut -f3 hg19_refGene.txt | sort | uniq -c | sort -k1,1nr | head -20 >tmp #建议存成文件，反复排序浪费CPU
awk '{print $1}' tmp | paste -s -d +
# -s：把文件的多行变成1行，默认分隔符是空格
# -d：改变分隔符
# awk 默认空格分隔，可有多个空格
## 7053+4938+4278+4256+4125+3780+3712+3670+3628+3465+3264+2902+2900+2855+2852+2639+2475+2299+1850+1562

awk '{print $1}' tmp | paste -s -d + | bc # bc：计算表达式

grep -w gene hg19_refGene.txt| wc 
awk '{if($3=="gene")print}' grep gene hg19_refGene.txt | wc #awk外边使用单引号，则内部需使用双引号

perl -alne '{print if $F[2] eq "gene"}' hg19_refGene.txt | wc

# 查看运行时间
time grep gene hg19_refGene.txt| wc 
time awk '{if($3=="gene")print}' hg19_refGene.txt | wc
time perl -alne '{print if $F[2] eq "gene"}' hg19_refGene.txt | wc

# 查看使用最频繁的命令
history | awk '{print $2}' | sort | uniq -c | awk '{print $1 "\t" $2}'|sort -k1,1n #注意uniq后的两列间的分隔符不规范，可能是多个不等的空格；awk后的单双引号问题
head * | less -NS #q

sort 常用参数总结

-k1,1：-k选项指定某列的排序方式。而每次使用 -k选项都要带上指定列的范围(start, end)。如果只指定一列，就为(start,start)了，像上面命令的 -k1,1就是。
-r: 逆向排序
-n: --numeric-sort 指定程序把列当做数值对待。
-g: 可按科学计数法进行排序
-V: --version-sort：在文本内进行自然版本排序，染色体按大小号显示
-t: --field-separator=分隔符使用指定的分隔符代替非空格到空格的转换